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Abstract. Second-order accurate upwind and 
centered schemes are presented in a framework that 
facilitates their analysis and comparison. The up- 
wind scheme employed consists of a reconstruc- 
tion step (MUSCL approach) followed by an up- 
wind step (Roe’s flux-difference splitting). The 
two centered schemes are of Lax-Friedrichs (L-F) 
type. They are the nonstaggered versions of the 
Nessyahu-Tadmor (N-T) scheme and the CE/SE 
method with e = 1/2. The upwind scheme is ex- 
tended to the case of two spatial dimensions (2D) 
in a straightforward manner. The N-T and CE/SE 
schemes are extended in a manner similar to the 
2D extensions of the CE/SE scheme by Wang and 
Chang for a triangular mesh and by Zhang, Yu, 
and Chang for a quadrilateral mesh. The slope 
estimates, however, are simplified. Fourier stabil- 
ity and accuracy analyses are carried out for these 
schemes for the standard ID and the 2D quadri- 
lateral mesh cases. In the nonstandard case of a 
triangular mesh, the triangles must be paired up 
when analyzing the upwind and N-T schemes. An 
observation resulting in an extended N-T scheme 
which is faster and uses only one third of the stor- 
age for flow data compared with the CE/SE method 
is presented. Numerical results are shown. Other 
improvements to the schemes are discussed. 

Introduction. When solving a fluid flow prob- 
lem, a researcher has the option of choosing between 
upwind and centered schemes using a quadrilat- 
eral or a triangular mesh. In this paper, trade-offs 
among these choices are discussed. Relations be- 
tween schemes and their strengths and weaknesses 
are shown, and improvements are suggested. 

The schemes employed are among the simplest 
second-order accurate schemes that can capture 
shocks and deal with unsteady problems. Both the 
upwind and centered schemes here use piecewise- 


linear reconstructions, i.e., MUSCL interpolants 
(monotone upwind schemes for conservation laws, 
Van Leer, 1977), which extend Godunov’s piece- 
wise constant method (1959). The key difference is 
that for the upwind scheme, numerical dissipation 
is added by the upwind step (Roe 1981, 1986), while 
for the centered schemes, dissipation is obtained by, 
loosely put, averaging the neighboring data (scheme 
ORD of Nessyahu and Tadmor 1990). 

The upwind step has a few drawbacks. Roe’s 
flux-difference splitting, which is mathematically 
rigorous and among the most popular, may cause 
oscillations as in the case of a slow-moving shock, 
or instability as in the carbuncle problem. The 
AUSM scheme (Liou and Steffen Jr. 1992, Wada 
and Liou 1997) does not have these problems, but 
it is not clear to this author why the scheme works. 
The upwind step is also sometimes perceived to be 
costly and difficult to grasp. In spite of these prob- 
lems, upwind schemes are popular because they 
work well for a large class of flows. The upwind 
step employed here is Roe’s splitting with an en- 
tropy correction described in Huynh (1995a). It 
can be derived by diagonalization and coded by 
stepping across one acoustic wave. The resulting 
scheme is concise and economical; the presentation 
below is also simpler than most presentations in the 
literature. Numerical solutions obtained with this 
upwind scheme for the ID Euler equations can be 
found in Huynh (1995a,b). The 2D extension of this 
scheme is conceptually straightforward. For other 
versions of upwind schemes, see, e.g., Barth and 
Jespersen (1989), Roe (1989), Hirsch (1990), and 
Venkat.akrishnan (1995). 

To avoid upwinding, a second-order accurate 
scheme, which extends the first-order scheme of 
Lax-Friedrichs (L-F), was introduced by Nessyahu 
and Tadmor (1990). There, the reconstruction step 
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is the same as that of the upwind scheme, but 
the upwind step is avoided by the use of a stag- 
gered mesh. The scheme employed here, however, 
is the nonstaggered version obtained by overlaying 
two staggered meshes to form a regular mesh. The 
drawback of the nonstaggered version is that the 
computing time doubles. The main reason this ver- 
sion is chosen is that the 2D extension retains its key 
advantage of simplicity especially when the mesh is 
triangular. In addition, the mesh and the boundary 
conditions can be chosen to be the same as those of 
the upwind scheme. The staggered version, on the 
other hand, requires two sets of meshes and two sets 
of boundary conditions — compromising the advan- 
tage of simplicity. Moreover, for a triangular mesh, 
the extension of the staggered version is quite in- 
volved (Arminjon et al. 1997) and, as will be dis- 
cussed, the advantage in computing time may no 
longer hold. Therefore, unless otherwise stated, we 
deal only with the the nonstaggered version below. 
Note that since there is no one-sided bias, the N-T 
scheme is centered. Also note that these L-F type 
centered schemes are different from the semidiscrete 
centered schemes made popular by Jameson et al. 
(1984). 

While the numerical dissipation of the N-T 
scheme is a lot less than that of the L-F method, 
it is still considerable. The CE/SE or conservation 
element and solution element method (Chang 1995) 
provides a way to adjust dissipation for these cen- 
tered schemes. Compared with the N-T scheme, 
the mesh, the balancing of fluxes, and the updates 
of the cell average quantities are essentially iden- 
tical. The difference is in the calculation of the 
slopes (of the linear interpolant). For CE/SE, the 
slopes must be stored, and due to the way the slopes 
are updated, numerical dissipation can be adjusted 
via a parameter called e. For the general CE/SE 
scheme (e / 1/2), the slope calculation is quite 
different from that in a typical MUSCL approach. 
When e = 0, the scheme has no numerical dissi- 
pation, i.e., it is reversible in time. Currently, the 
CE/SE member employed in essentially all practi- 
cal calculations corresponds to e = 1/2. For this 
reason, we restrict our attention to this member 
and, from this point on, unless otherwise stated, 
the term CE/SE refers to the member with e = 1/2. 
Note that there are numerous differences in termi- 
nology between (Nessyahu and Tadmor 1990) and 
(Chang 1995); here, the terminology in the former 
is often employed. 

The CE/SE schemes were extended to 2D for un- 
structured triangular meshes by Wang and Chang 


(1999) using the nonstaggered version. What is 
novel about this extension is that the spatial do- 
main where each reconstruction is valid at the be- 
ginning of the time level is a hexagon, which is 
roughly twice as big as the triangular cell. This 
CE/SE approach to extension is also applied here 
to the N-T scheme. (Such an extension was men- 
tioned as a coupled version for the CE/SE scheme 
in Chang et al. (1999) and has recently been 
incorporated — independently from this work — as 
an option in the CE/SE code; private commu- 
nication with Drs. Ananda Himansu, Ching Y. 
Loh, and Xiao- Yen Wang. Note, however, that 
the extended N-T scheme presented here has nu- 
merous differences resulting in a scheme which 
is faster and requires considerably less storage.) 
The quadrilateral-mesh extension for the CE/SE 
method can be found in Zhang et al. (2002); see 
also Cook (1999). The CE/SE schemes have been 
applied to solve numerous practical problems in two 
and three dimensions with a lot of success, espe- 
cially in aeroacoustics. 

For a structured quadrilateral mesh, in a manner 
similar to the ID case, the nonstaggered mesh in 
Zhang et al. (2002) can be obtained by overlaying 
two staggered meshes in Arminjon et al. (1995) 
and Jiang and Tadmor (1998) (see also Jiang et al. 
1998). For a triangular mesh, however, a similar 
statement does not hold; in fact, the nonstaggered 
extension in Wang and Chang (1999) appears to 
have numerous advantages over a staggered-mesh 
extension (remark (c) in §6 below). 

In this paper, the schemes involved are first pre- 
sented for the ID advection equation where key 
ideas and trade-offs can already be seen. It is shown 
that the N-T and the CE/SE (e = 1/2) schemes are 
respectively the centered counterparts of Van Leer’s 
first and second upwind schemes. Next, the exten- 
sions of the upwind, N-T, and CE/SE schemes to 
the ID Euler equations are explained. Then, exten- 
sions to the 2D Euler equations on a quadrilateral 
and a triangular mesh as well as the simplified slope 
estimates are described. Comparison of schemes via 
Fourier stability and accuracy analyses are carried 
out. Here, for a triangular mesh, we must pair up 
the downward and upward pointing triangles when 
analyzing the upwind and N-T schemes. 

Concerning the two L-F type methods, the N-T 
and CE/SE schemes are shown to produce essen- 
tially the same numerical solutions. The former has 
the advantage of better coupling; consequently, it 
converges better for a steady state problem. In ad- 
dition, the following observation yields an extended 
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N-T scheme which is faster and requires consider- 
ably less storage than the CE/SE method (slopes 
need not be stored): instead of gathering fluxes to 
update the solution, fluxes are distributed to the 
neighboring cells. This observation — made possible 
in part by the simplified slope estimates — results 
in each cell needing to be visited only once during 
each time step instead of three or four times as in 
the case of the CE/SE scheme (three for a triangu- 
lar mesh, and four for a quadrilateral mesh). 

This paper is essentially self-contained. Readers 
who are interested only in the Euler equations can 
start with §4. The paper is organized as follows. 
In §1, the first-order accurate upwind and centered 
schemes are presented as two different ways of sta- 
bilizing the Euler forward scheme. In §2, upwind 
schemes are developed for the ID advection equa- 
tion using piecewise linear reconstructions. The 
centered counterparts of these upwind schemes are 
described in §3; here, a comparison via Fourier anal- 
ysis is carried out. Section 4 deals with extensions 
to the case of the ID Euler equations; §5, the 2D 
Euler equations on quadrilateral meshes; and §6, 
on triangular meshes. In §7, Fourier analysis for 
the 2D case is carried out. Numerical examples are 
shown in §8, and conclusions are presented in §9. 

The author wishes to thank Prof. B. van Leer, 
Mr. C. Steffen, Jr., Drs. S.-C. Chang, R. Chima, A. 
Himansu, D. Jacqmin, P. Jorgenson, M.-S. Liou, C. 
Y. Loh, A. Suresh, and X.-Y. Wang for several illu- 
minating discussions. This work was supported by 
the Engine System Noise Reduction project, which 
is part of the Quiet Aircraft Technology Program 
at the NASA Glenn Research Center. 


where t is time, x distance, and uq(x) the initial 
condition. By assuming that uo(x ) is periodic, 
boundary conditions are straightforward and are 
omitted. The exact solution is 

u(x, t) = uq(x — at). (1.2) 

To discretize the above problem, let h be the 
mesh spacing, and Xj = jh, j = 0, 1, 2, ... be a uni- 
form mesh. Let At be the time step and t n = nAt 
be the time level. At time t n , let w” be an approx- 
imation to the the solution u at Xj. Assume that 
we know u” for all j ; we wish to calculate u" +1 . 
To simplify the notation, the superscript n in w" 
is omitted and the data is denoted by Uj\ all other 
superscripts, however, are retained. 

Next, set 

a = aAt/h. (1.3a) 

Then, the quantity |oj is the Courant number. As- 
sume that the time step satisfies the CFL (Courant, 
Friedrichs, and Lewy) condition 

\a\ < 1; (1.3 b) 

loosely put, information propagates no more than 

one mesh size per time step. 

1.1. Euler forward scheme. This scheme is 
given by the FTCS differencing, 

Uj +1 =uj -<j\{u j+ i -Uj-i). (1.4) 

Again, u ” is abbreviated to Uj. 

For the Fourier (or Von Neumann) stability anal- 
ysis, set Xj — j and 


1. First-order accurate schemes for advec- 
tion. The first-order case is straightforward, but it 
conveys the ideas and the trade-offs. In addition, 
improvements and nonstandard observations will be 
made. The simplest discretization — forward-time 
centered-space (FTCS) — leads to the Euler forward 
scheme, which is unstable. To stabilize this scheme, 
we need numerical dissipation. Dissipation can be 
added by using an upwind-biased difference, which 
results in an upwind scheme, or by using a dissi- 
pative start-off value, which results in a centered 
scheme. 

For simplicity, consider the advection equation 
with constant speed a, 


du du 

m +a d^ = 0 ' 

u(x,0) = u 0 {x), 


(1.1a) 

( 1 . 16 ) 


where I = \/— I and w is the wave number, — tt < 
w < 7T. Equation (1.4) implies 

u n+i = e ijw y l _ (J 1 ( e Iw _ e -iw ) ] 

The quantity in the square brackets is the amplifi- 
cation factor of the Euler forward scheme: 

A = l-<r[\{e Iw -e~ Iw )\. (1.5) 

If a > 0, then, the CFL condition (1.3b) implies 
0 < a < 1. Next, the slop e u x by central difference 
at j = 0 is given by \ C& = \(e Iw — e~ Iw ) shown in 
Fig. 1.1, where the point C corresponds to j = — 1, 
A to j =0, and E to j = 1. Thus, the amplification 
factor lies on the line segment AB, where A corre- 
sponds to a = 0 and B, a = 1. Since this range is 
outside the unit circle, the scheme is unstable. Also 
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Fig. 1.1. Fourier Analysis. The circle is of radius 1 
in the complex plane. The amplification factor for 
the Euler forward scheme lies on the line segment 
AB; that for the first-order upwind scheme, AC; 
and that for the L-F scheme, DC. To reduce nu- 
merical dissipation for the L-F scheme, instead of 
D, we can employ a start-off value which blends D 
and A, e.g., (1 — \<r\)uj + \a\ +Uj+i). 

note that after a time step corresponding to o , the 
exact solution at j = 0 is e~ Iaw . The amplification 
factor is an approximation of this exact solution. 

1.2. First-order accurate upwind scheme. 

The Euler forward scheme can be stabilized by an 
upwind-biased difference for ( u x )j . The result is 
(Courant, Isaacson, and Rees, 1952), 

u n +i _ f u i ~ &{uj — uj- 1 ) if a > 0, /j g, 

l \uj - a(u j+1 -Uj) otherwise. ' ' 

The amplification factor for the first-order upwind 
scheme is 

(l-a(l-e- Iw ) if a > 0, ( , 

\ 1 — a(e Iw — 1) otherwise. 

If a > 0, then 0 < a < 1, and the amplification 
factor lies on the line segment AC shown in Fig. 1.1. 
Since this range is inside the unit circle, the scheme 
is stable. 

Notice that if the wind direction is incorrectly 
determined — for systems of nonlinear equations, we 
may potentially run into this problem due to vari- 
ous approximations — the solution would follow the 
direction EA but would lie outside the circle and, as 


a result, the scheme may encounter stability prob- 
lems. 

1.3. First-order centered scheme (L-F). 

Another way to stabilize the Euler forward scheme 
is by employing a dissipative start-off value: replac- 
ing Uj in (1.4) by the average of the two neighbor- 
ing values. The result is the Lax-Friedrichs (L-F) 
scheme (Lax 1954): 

u j +1 = H u J-i + u j+i) ~ a U u i + 1 - C 1 - 8 ) 

Note that when At = 0, the solution, instead of 
being uj, is |(uj_i 4- «j+i). Such a quantity is 
called a start-off value here. 

The amplification factor for the L-F scheme is 

A = \{e Iw + e~ Iw ) - o ±(e Iw - e~ Iw ). (1.9) 

If a > 0, then 0 < a < 1, and the amplification 
factor lies on the line segment DC shown in Fig. 1.1. 
As a result, the scheme is stable. Observe that the 
two key advantages of the L-F scheme are simplicity 
and stability (the solution lies well inside the unit 
circle) . 

Concerning the disadvantages, by (1.8), the so- 
lution at j does not depend on the data Uj , but 
depends only on the data at the two neighboring 
indices. Consequently, the scheme has the odd-even 
decoupling problem. This problem can also be seen 
via the staggered-mesh formulation. Indeed, sup- 
pose at time level n, we have data at even indices 
only. At time level n + 1, (1.8) yields the solu- 
tion at all odd indices. At time level n + 2, we 
can obtain the solution at all even indices again, 
etc. Such a scheme employs a staggered mesh. The 
other staggered-mesh solution has data at odd in- 
dices at time level n. Thus, the regular-mesh ver- 
sion (1.8) above consists of two completely indepen- 
dent staggered-mesh solutions. As a consequence, 
the L-F scheme needs twice the number of mesh 
points to have the same resolution as the first-order 
upwind scheme. (This observation can also be seen 
from Fig. 1.1.) In other words, the convenience of a 
regular mesh is achieved at the cost of a) twice the 
amount of calculations and b) the odd-even decou- 
pling problem. 

The next two observations hold for the L-F as 
well as the second-order L-F type schemes in §3. 
First, in an unsteady calculation, if we wish to an- 
imate the solution, we should use the solutions at 
even time steps only (or odd ones only). Similarly, 
to check for convergence to a steady state, one must 
compare the solution at time level n with that at 
time level n — 2. Second, when o = 0, the L-F 
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scheme has no phase error while the damping er- 
ror reaches a maximum. Consequently, in practical 
calculations, it is desirable to employ a time step as 
large as possible to minimize damping. 

1.4. Scheme blending Lax- Friedrichs and 
Euler forward. The two drawbacks of the L-F 
scheme (odd-even decoupling and too much numer- 
ical dissipation) can be dealt with by blending it 
with the Euler forward method. In other words, in 
Fig. 1.1, instead of starting off at the point D, we 
can employ a start-off value which blends D and 
A: the quantity i + Uj + i ) in (1.8) is replaced 
by (1 — |<r|)ttj + | cr | |( Uj-i +Uj+i). The resulting 
scheme is 

u” +1 = [(1 - \o\)uj + | cr | 4-tij+i)] 

- \ (U j+ i -Uj- 1). 

( 1 . 10 ) 

The above centered scheme turns out to be iden- 
tical to the first-order upwind scheme. For the 
case of second-order accuracy, however, the blended 
scheme is different from the second-order upwind 
one. 

Note that in extending (1.10) to the Euler equa- 
tions, the quantity \a\ takes the form \ A\At/h where 
A is the Jacobian matrix. For simplicity, |A| can be 
replaced by the sum of the magnitude of the local 
flow speed and the speed of sound. 


2. Upwind schemes via MUSCL approach. 

The above schemes were derived by finite differenc- 
ing. To ensure that shocks are captured, fluxes need 
to be balanced via the integral form of the equation 
(Lax 1954). Oscillations near shocks can be sup- 
pressed by using linear functions (or polynomials) 
to approximate the data in each cell and by limiting 
the slopes of these linear functions so that they are 
not too steep near a discontinuity (Van Leer 1977). 

Integrating (1.1a) on the interval [a,/?], one ob- 
tains 


d_ 

dt 



dx + au(f3 , t) — au(a , t) = 0, 


( 2 . 1 ) 


where au(a, t) is the flux (or rate of flow) of u across 
interface a at time t. 

Next, let Xj+ 1/2 = (j + 1/2 )h be the cell inter- 
faces and Xj = jh the cell centers of a uniform 
mesh, j =0,1,2,... For each cell [a;j_i/ 2 ,ar J+1 / 2 ], 
second-order accuracy is obtained by applying the 
midpoint rule to (2.1): 


i + ¥<<“£$ -<“£$); ( 2 - 2 > 



Fig. 2.1. Interface fluxes for second-order accurate 
upwind schemes. At each interface j + 1/2 time 
level n + 1/2, the reconstruction rj to the left of 
the interface yields u L while rj + 1 yields u R . The 
upwind value is employed for the flux calculation. 


here, uj (i.e. , u’f ) approximates the average value of 
u in the cell j at time t n , and \ approximates 
the value at the interface j -4-1/2 at time level n + 
1/2. Assume that we know the data Uj (i.e., w") 
for all j ; we wish to calculate the solution «" +1 . 

2.1. First-order upwind scheme. For the 

first-order case, in each cell j, the data are approx- 
imated by a constant function: rj(x) = Uj. Next, 
applying the FTCS (forward-time centered-space) 
approximation to (2.1), one obtains 

u] +1 = uj + ^{auj_ 1/2 - au j+1/2 ). (2.3) 

At each interface j + 1/2, the data to the left is the 
constant function Uj , that to the right, Uj + 1 . The 
flux is given by the upwind choice 


au j+ 1/2 — 


auj if a > 0, 
auj+i otherwise. 


(2.4) 


For consistency with the case of the Euler equations 
later, at each interface, a and u are combined as 
above. 


2.2. Second-order upwind scheme. For 

second-order accuracy, the data are approximated 
by a linear function in each cell. To calculate the 
flux in (2.2), suppose {(u x )j} are known. 

They can be approximated by, e.g., the central dif- 
ference (also called the average slope) 

(' u x)j = 2 h( u j+l — u j-l)- (2-5) 

(A weighted average for (u x )j will be discussed 
later.) The time derivative ( u t )j follows from the 
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advection equation, 


(' u t)j — a ( u x)j- (2-6) 

In the domain [xj_ i/ 2 ,^j+i/ 2 ] x [i T \i n+1 ], the solu- 
tion can then be approximated by a linear function 
rj called the reconstruction , 

Tj{x. t) = Uj + (u x )j(x — xj) + (u t )j(t — t ). (2.7) 


Note that the term ‘reconstruction’ is typically used 
for rj(x,t n ). While rj (x. t ) was employed by Harten 
et al. (1987) and Nessyahu and Tadmor (1990), 
it was not given a name. Here, we use the term 
‘reconstruction’ in the extended sense above. 

At the two interfaces of cell j , we can march to 
the half time step using rj . At each interface Xj + 1 / 2 
and time t n+1/l2 , we now have two values. The value 
from cell j is denoted by u L ; that from j + 1, by 
u R , as shown in Fig. 2.1, 

u l =u j + | ( u xh + 

U R = U j + 1 — 2 ( U x)j+1 3 2~{ U t )j+l- 

The upwind choice for the flux is 


n+ 1/2 J au L if a > 0, 
au j+ 1/2 | aUR otherwise. 


( 2 . 8 ) 


This flux and (2.2) complete the description of the 
second-order upwind scheme. 

The calculation (2.6) of the time partial deriva- 
tive from the spatial one follows Hancock’s obser- 
vation (Van Albada et al. 1982). It simplifies the 
extension to systems of equations. Also notice that 
the Taylor series rj(x,t) yields a result identical to 
that by the method of characteristics (1.2) applied 
to the linear initial condition rj(x,t n ). 

Note that if a = 0, then there is ambiguity in 
defining the interface value in (2.8), but 

since the flux is zero, this ambiguity causes no prob- 
lem. In fact, the flux can be expressed without a 
conditional statement: 


au "+i /2 = U au L + au R ) - ±\a\(u R - u L ). (2.9) 

Similarly, for systems of equations, there is no ambi- 
guity in splitting the flux, and the technique works 
well. Splitting the variables, however, may result 
in unwanted oscillations. Also note that the term 
— \\ a \{u R — u L ) produces the upwind-biased effect 
and thus numerical dissipation. 

If a > 0, expressions (2.2) and (2. 5-2. 8) yield the 
following solution 

= Uj + <r(wj_i — uj) + 

\cr{\ — rr)[uj + Uj-\ — Uj - 2 — «j+ i]- 


This scheme was formulated via finite differencing 
by Fromm (1968). The above reconstruction formu- 
lation follows that of (Van Leer 1977) except that 
the method of characteristics there is replaced by 
the Taylor series (2.7) here. The scheme is the first 
and also the least accurate among a series of five 
schemes in the cited reference. 

Fourier analysis yields the following amplification 
factor for the second-order upwind scheme, 


A - 1 + a(e~ Iw - 1) + 

j<t(1 — <t)(1 + e~ Iw — e~ 2Iw — e Iw ). 

(2.10) 

To avoid oscillations near a discontinuity, instead 
of the average (2.5), a weighted average (Van Al- 
bada et al. 1982) can be employed for the slope. 
For any two real variables x and y, define 


wtav ( x , y) 


x 2 y -I- y 2 x 
x 2 + y 2 + 10 -20 


( 2 . 11 ) 


If x/y k, 1, the above expression yields essentially 
the average \(x + y); on the other hand, if x/y ss 0 
or oo, the result is essentially the value with smaller 
modulus. The slope formula takes the form 


{u x )j = £ wtav(ttj + i - Uj , Uj - Uj- 1 ). (2.12) 

Note that between two slopes \{uj — ttj-i) and 
\(uj + 1 — Uj), the above weighted average produces 
a result which is biased toward the less steep slope. 
Since such a result is closer to the zero slope of the 
first-order upwind scheme than the result by the av- 
erage (2.5), the above weighted average adds a con- 
siderable amount of numerical dissipation. In fact, 
it yields a scheme which is only first-order accurate 
near an extremum — a well-known drawback of all 
standard weighted averages or limiter functions. 

2.3. Scheme II of Van Leer. Next, we de- 
scribe an upwind scheme which carries along (stores 
and updates) the interface values. The scheme is 
the second one presented in (Van Leer 1977) ex- 
cept, again, the method of characteristics employed 
there is replaced by the Taylor series (2.7) here. 
While this upwind scheme does not extend easily 
to the case of the Euler equations (more on this 
later), its centered counterpart, which turns out to 
be the CE/SE scheme with e = 1/2, does. 

Suppose the data consist of the cell averages {uj} 
as well as the interface values {uj +1 / 2 }- We wish to 
calculate {u/ +1 } and First, instead of the 

average (2.5), the slope is defined via the interface 
values: 


( u x)j — f l ( u j+i/2 u j- 1 / 2 )- (2-13) 
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With the above slope formula, the cell averages 
{u" +1 } are calculated as before, by (2.2) and (2.6)- 
(2.8). The interface values are updated by using 
the reconstruction in the upwind cell: 


n -\- 1 

3 + 1/2 ~ 


f rj(x j+1/2 ,t n+1 ) 
\r j+1 (x j+1/2 ,t n+1 ) 


if a > 0, 
otherwise. 


(2.14) 


This completes the description of the scheme. 

The above scheme has the same leading (phase) 
error compared with the scheme using (2.5). The 
dissipation error, however, is reduced by two thirds. 
Thus, carrying along the interface values improves 
the dissipation error considerably. 

Note that extending (2.14) to the Euler equations 
appears to be difficult: when a = 0, there is no 
ambiguity in defining the flux, but the choice of the 
interface value becomes ambiguous. 


3. Centered schemes via MUSCL appoach. 

Each upwind scheme in §2 has a centered counter- 
part in this section. For centered schemes, the spa- 
tial region where the reconstruction is valid at time 
t ~ t n . which is called the reconstruction cell here, 
is twice the size of the original cell. Together, these 
reconstruction cells cover the domain twice whereas 
the original cells cover the domain once only. 

More precisely, on the same mesh of cells 
^b+1/2] as i n the upwind case, let the recon- 
struction cell j be the interval [xj_i,Xj + i] of length 
2h defined by the centroids of the two neighboring 
cells. Assume that we know the data Uj (i.e., u ”) 
for all j; Uj approximates the average value of u on 
the reconstruction cell j. We wish to calculate the 
solution u" +1 . See Fig. 3.1. 

3.1. First-order centered scheme (L-F). 

Here, we rederive the L-F scheme (1.8) from the 
MUSCL perspective. In each reconstruction cell j, 
the data are approximated by a constant function: 
rj(x) = uj. Next, when updating the solution at j, 
the reconstruction rj is ignored , while rj- 1 and rj + 1 
are employed. See Fig. 3.1. Thus, the start-off (av- 
erage) value is 

u *j = +«j+i)- (3- 1 ) 

Applying the FTCS approximation to (2.1) on the 
reconstruction cell j, one obtains 

Uj + = 2 {uj-l +Mj+l) + 2^{ au j-l — au j+l)- (3-2) 

Here, there is no need of upwinding because the flux 
auj - 1 is evaluated at the center of the reconstruc- 
tion cell j — 1, and the flux auj + 1, the center of the 
reconstruction cell j + 1. 



Fig. 3.1. Centered schemes. The interval 
[xj-\,Xj + i\ is called the reconstruction cell j (the 
spatial region where the reconstruction is valid at 
time t = t n ). The two adjacent reconstruction cells 
j and j + 1 overlap on [xj,xj+ 1]. For each j, 
the reconstruction rj(x,t) is valid on the triangle 
defined by the three corners (, xj-i,t n ), (xj+i,t n ), 
(xj,t n+1 ). When calculating uj +1 , the reconstruc- 
tion rj is ignored, while r,-_i and fj+i are em- 
ployed. For second-order accuracy, the values at 
the points marked by x and * are needed. They 
are uj - 1 + ^(«<)j- 1, u j+ 1 + 1, and 

Uj — 1 T ^(.Ux'lj— 1? and Uj + 1 ^(ttaOj+i* 

The next remark concerns the initial values for 
the L-F scheme. The remark also holds for the N-T 
and the CE/SE schemes below. Suppose the ini- 
tial condition is a step function, say, Uq(x) = 1 for 
x < 0, and uq(x) = 0 otherwise. Suppose the cell 
centers are Xj = j, j = -N, ...,7V where TV = 50. 
Then the solution at a later time t n in standard 
textbooks exhibits a staircase pattern: 14 = u 2 , 
and us = 14, etc. One way to avoid this problem is 
to define the initial data as the average of uo(x) on 
the reconstruction cell j, i.e., = 1/2. 

3.2. Second-order centered scheme (N-T). 

This scheme is the centered counterpart of the 
second-order accurate upwind scheme of §2.2. It 
was presented by Nessyahu and Tadmor (1990) us- 
ing limiter functions of minmod type to define the 
slope and was named ORD. For consistency with 
the upwind scheme here, Van Albada’s weighted 
average is employed. In addition, the average slope 
(2.5) makes Fourier analysis possible. 

As in the upwind case, the slopes are given by 


(u x )j — 2 h (llj+i Uj- 1), 

(3.3) 

i U t)j = ~ a ( u x)j- 

(3.4) 


At time t n , on the reconstruction cell j, the data is 
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assumed to be the linear function Uj+(u x )j(x — xj). 
Since information propagates a distance no more 
than h per time step, information from outside the 
reconstruction cell has not reached its center regard- 
less of the wind direction. As a result, the above ut 
is valid in the triangle defined by the three corners 
(xj_i,f ra ), ( Xj+i,t n ), (xj,t n+1 ). In this region, the 
solution can be approximated by, as in (2.7), 

rj(x,t) =Uj + ( u x )j{x - Xj) + (u t )j(t -t n ). (3.5) 

When updating the solution at j. again, the re- 
construction rj is ignored, while rj-i and rj + 1 are 
employed. The start-off (average) value is 

u j = 2 [ U J- 1 "t" ^( W ®);7-1 J r' u o+ 1 — 2 ( w z)j+l]- (3-6) 

See Fig. 3.1. The midpoint rule for (2.1) on the 
reconstruction cell j yields 

u T l = u j + f ( au J-i /2 - «<+i 1/2 )- ( 3 - 7 ) 

Note that the use of Uj instead of uj results in 
numerical dissipation which stabilizes the scheme. 
The values and 2 are evaluated by the 

reconstruction function: for any j, 

u j+l/2 =Uj + At ( Ut )j . (3.8) 

The above completes the description of the N-T 
scheme. 

With expression (3.3) for the slope, the result is, 

U ] +1 = + \(Uj - Uj-i) + 

— 4 ( u j+2 — Uj) ] + (3 9) 

\ ° { [Uj - 1 - \a{uj - Uj- 2 )] - 

[ u j + 1 — j (J ( u j +2 ~ u j )] }• 

At first glance, since the right hand side above 
involves Uj, the N-T scheme appears to be (odd- 
even) coupled. This coupling, however, is not a 
strong coupling, as can be seen in the following ex- 
ample. Suppose the data at odd indices are 1 and 
those at even, —1. Then the solution by the N-T 
scheme after two time steps remains identical to the 
initial data for all time step sizes. It is desirable for 
a scheme to damp out this odd-even data. In spite 
of this observation, based on the author’s experi- 
ence, the N-T scheme produces solutions with no 
odd-even noise in practice provided that the ini- 
tial condition is appropriate as discussed in the last 
paragraph of §3.1. 

The amplification factor of the N-T scheme is 
given by, after some algebra, 

A = cos(w) — /<jsin(w) + | (1 — a 2 ) sin 2 (w ). (3.10) 


The next remark concerns the staggered-mesh 
version. At time level n, suppose we have data at 
even indices only, and we wish to calculate the so- 
lution at odd indices. Then, the slope at an even 
index 2 j is (« 2 j +2 — t< 2 j- 2 ) / (4/i). For the non- 
staggered version above, the slope at index 2 j is 
(w 2 j+i — U2j-i)/(2h). Because the latter slope uses 
values at locations closer to j , the numerical dissi- 
pation of the nonstaggered version is only 1/3 that 
of the staggered version. The phase errors, however, 
are essentially the same (for small wave numbers). 
Thus, the nonstaggered version costs twice as much 
as the staggered one, but its numerical dissipation 
improves by two thirds. 

As in the first-order case, to reduce numerical 
dissipation for the nonstaggered N-T scheme when 
<7 is small, we can use a blended start-off value: with 
u* defined by (3.6) 


U" +1 = {(1 - W\)Uj + \a\u*} + 

1 a 4 -( n-t~l/2 n-\-l/2\ ^ ^ 

t(auj_/ - au j+1 ). 

The resulting scheme also damps the data of 1 at 
odd and -1 at even indices. 

3.3. CE/SE scheme. The centered counter- 
part of Van Leer’s second scheme in §2.3 turns out 
to be Chang’s CE/SE scheme (the e = 1/2 mem- 
ber). To show this fact, we describe the CE/SE 
scheme as one which carries along (stores and up- 
dates) the slopes, and then as one which carries 
along the interface values. 

At time £”, suppose the data {uj} as well as the 
slopes |(wz)j} are known and stored. We wish to 
calculate {w" +1 } and {(uz)” +1 }- The cell average 
quantities {?.<" 1 1 } are updated exactly as in the N-T 
case: by (3.4)-(3.8). To update the slope, first, for 
each j , the point value (as opposed to cell average) 
m" + 1 is calculated by the reconstruction: 

u" +1 = rj(xj,t n+1 ) — Uj + (u t )jAt. (3.12) 

The slope is updated by the two neighboring point 
values: 


(u*)" +1 = Jh^j+i (3-13) 

If the weighted average is employed, the cell average 
u] +1 is also needed: 

« = s wtav(ff”+ 1 1 - u]+\u] +1 - u n j + l). 

Clearly, instead of storing the slopes {{u x )™ +1 }, 
we could store the interface values {u" +1 }. Such a 
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scheme takes the following form: at time t n , sup- 
pose the data { Uj } as well as the point values { uj } 
are known and stored. To update the solution, first, 
the slope is given by 

(u x )j = £ wtav(itj +1 - Uj ,Uj — Sj-i). (3.14) 

The cell average u' ,+ 1 is then updated by again 
(3.4)-(3.8). The point value fi” +1 is updated by 
(3.12). 

Note that the CE/SE scheme is odd-even decou- 
pled, as can be seen by the following staggered-mesh 
formulation. At time level n, suppose we have data 
U2j at even indices and U2j+i at odd indices only. 
Here, {u2j+i} are the values at the interfaces of the 
cells [x 2 j~i, X 2 j+i\- At time level n + 1, we calculate 
the solution u^j+i an< l 1 ■ The other staggered- 
mesh solution is obvious. Thus, the nonstaggered 
version is a result of overlaying two independent 
staggered-mesh solutions. 

The above observation also shows that for ID 
advection, in the form that carries along the inter- 
face point values, the CE/SE scheme is the centered 
counterpart of Van Leer’s second scheme described 
in §2.3. 

Notice that the CE/SE scheme in above form 
needs less storage: for the 3D case, instead of stor- 
ing u, u x , u y , and u z , we only need to store u and 
(the point value) u. If we wish to adjust numerical 
dissipation (e ^ 1/2), however, we have to store the 
slopes. 

A few remarks concerning the CE/SE and the 
N-T schemes are in order. First, if the slopes are 
carried along, the CE/SE scheme has a very com- 
pact stencil: when updating the solution at cell j, 
only the data at the immediate neighbors j — 1 and 
j + 1 are employed. But the slopes at j — 1 and j + 1 
use the point values at j — 2 and j + 2. Thus, in 
the form of carrying along the interface point val- 
ues, the stencil of the CE/SE scheme is the same as 
that of the N-T method. 

Next, after one time step, the cell average update 
u” +1 has an error of 0(h 3 ) (the solution is exact if 
the data are on a parabola). The point value update 
u n + x , by (3.12), has an error of 0(h 2 ) (exact when 
the data are linear). Therefore, when calculating 
the slopes (u x )™ +1 , the cell averages u"** and u"^ 1 
appear to be a better choice than the points values 
u’J+l and Uj+i- It turns out, however, that the 
point values lead to a scheme with the same phase 
error but less dissipation error (more details later), 
i.e., the CE/SE scheme (e = 1/2) is less dissipative 
than the N-T scheme. 


For Fourier stability analysis, set 



Then, the solution vector of the CE/SE scheme can 
be written as 


u”« = 

Bu^i + Cu j+1 , 

(3.17) 

where 



B = + 

1(1 — o-) (1 + er)^ , 

(3.18) 


and 

C = (^ 1 r a) -i( i -jo (1+ff) ). (3.i9) 

The two eigenvalues of the matrix e~ ,w 'B+e Iw C 
are the amplification factors of the CE/SE scheme. 
After some algebra, 

A^ = I cos(w) — 1(7 sin(w) ± 

, (3.20) 

| \J 1 + (1 - 2a 2 ) sin 2 (w). 

Here A + , which approximates e~ Iaw and deter- 
mines the accuracy of the scheme, is the principal 
amplification factor, while A~ is the spurious one. 

Our next remarks concern the start-off slopes and 
the reversible CE/SE scheme. This scheme involves 
looking both forward and backward in time (Chang 
1995). It is described below using a start-off slope 
and a forward only time evolution. For the case 
of the Euler equations, such a forward marching 
scheme provides an approximation to the reversible 
scheme. 

At time t n , suppose {uj} and {(u^-} are known. 
Then if At = 0, the slope (u x )™ +1 of the above 
CE/SE scheme (e = 1/2) reduces to the start-off 
slope 

= (3.21) 

For At > 0, the slope update (again, e = 1/2) is 
given by 

K)" +1 = K)* + f [KW i - Mi- 1]. (3-22) 

The reversible CE/SE scheme (e = 0) can be de- 
scribed as a scheme with a different start-off slope. 
Instead of (u x )*, the start-off slope is defined by the 
values employed in the start-off average (3.6): 

[iP‘ X )j] R = {uj+i — 2 i u x)j+l ~ [ u j-i+ 2 ( u z)j-i]}- 

(3.23) 
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Fig. 3.2. Dissipation errors (per time step) as func- 
tions of a for w = 7r/12 (relatively low frequency). 
For each scheme, this error is proportional to w 4 . 
When (7 = 0, the N-T and CE/SE schemes have 
a maximum amount of numerical dissipation. For 
reference, the error of the Lax- Wendroff scheme is 
also plotted. While the Lax- Wendroff scheme has 
relatively small errors, it causes oscillations near 
shocks. 

The slope update for the reversible scheme is 
Oz)" +1 = [(«*),•]]? + 1 - (u t )j- 1]. (3.24) 

To obtain the general CE/SE schemes, we simply 
blend the two start-off slopes (3.21) and (3.23) us- 
ing a parameter e by [(«x)j]R + 2e{(u x );f -[(Ux)j]r}. 
The evolution for the slope, i.e., the quantity 
j£[(ut)j+i — («t)j_i], is the same for all e. As 
such, the CE/SE schemes have an adjustable start- 
off slope. 

It is remarkable that the start-off average u* via 
(3.6) produces numerical dissipation, whereas the 
start-off slope can have the opposite effect: it can 
reduce dissipation. Note that one can derive the 
upwind counterparts of the reversible and the gen- 
eral CE/SE schemes (0 < e < 1) as well. These 
schemes, however, are beyond the scope of this pa- 
per. 

Also note that the reconstruction cell here is iden- 
tical to the spatial part (or spatial projection) of the 
solution element in Chang (1995), and the control 
volume on which fluxes are balanced here, the con- 
servation element there. 

3.4. Accuracy comparison. The upwind and 

centered schemes described above are stable for 
| <7 1 < 1. To compare accuracy, recall that the am- 
plification factor A approximates the exact amplifi- 
cation factor e~ Iaw . The phase error per time step 
is Arg(A) — (—aw). For the second-order schemes 
discussed here, this quantity is proportional to w 3 
when the wave number w is small. The dissipation 



Fig. 3.3. Phase errors (per time step) as functions 
of a for w = 7t/12. For each scheme, this error is 
proportional to w 3 and is the leading error. Note 
that the phase errors of the two centered schemes 
are nearly identical. (They are identical on a simi- 
lar plot for w = 7t/32./ 

error per time step is \A\ — 1, which is proportional 
to w 4 . 

Figure 3.2 shows the dissipation errors as func- 
tions of (7 for w = 7t/ 12 (relatively smooth data). 
For ease of reference, the error of the Lax- Wendroff 
scheme is also plotted. 

Figure 3.3 shows phase errors as functions of a for 
again, w = 7t/12. The phase errors of the centered 
schemes are nearly identical. Notice that when a = 
1/2, the upwind scheme has no phase error, while 
its dissipation error reaches a maximum. 

Observe that for ID advection, the upwind 
scheme is considerable more accurate than the non- 
staggered N-T and CE/SE schemes. The reason 
is that the cells of the upwind scheme are half the 
size of the reconstruction cells of the two centered 
schemes. 

Between the N-T and the CE/SE schemes, 
Fig. 3.2 shows that the latter has less numerical dis- 
sipation. For a small CFL number, however, both 
schemes have about the same amount of dissipation. 
In practice, the weighted average adds additional 
dissipation and the results by these two schemes 
are nearly the same, as will be shown. 

4. Schemes for the Euler equations. The 

second-order upwind, N-T, and CE/SE schemes are 
extended to the case of the Euler equations in this 
section. For this case, the N-T and CE/SE schemes 
remain simple; in fact, the equations are identical to 
those in the case of advection except that the sym- 
bols are boldfaced. As for the upwind scheme, we 
need additional techniques, namely, Roe’s splitting 
with an entropy correction. Roe’s method consists 
of a diagonalization and an upwind side selection 
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in the characteristic frame. The entropy correc- 
tion serves the purpose of excluding non-physical 
solutions and avoiding kinks (glitches) in the solu- 
tion near sonic points. The entropy correction in 
Huynh (1995a), which amounts to adding dissipa- 
tion when the wave speeds spread past zero, is em- 
ployed. While these concepts for the upwind step 
are somewhat involved, the coding remains concise 
and economical. The derivation for this upwind 
step is carried out below to show the complexity — 
or simplicity — of upwinding. Also note that the 
presentation of Roe’s scheme here is simpler than 
most presentations in the literature. 

The one-dimensional flow of an inviscid and com- 
pressible gas obeys the conservation laws for mass, 
momentum, and energy: 

U t + F(U) x =0, (4.1) 


U = ( A ] , F = lpS+pY 

\ e J \ (e+p)uj 

where t is time, x distance, p density, e total energy 
per unit volume, u velocity, and p pressure. Let 7 
be the ratio of specific heats. Then for a perfect 
gas, 

P = (7 - l)(e - \pu 2 ). (4.2) 

Integrating (4.1) on the interval [a,/?], one ob- 
tains the integral form 


d_ 

dt 



<te + F(U(0,t)) -F(U(a,<)) 


= 0, 


(4.3) 

where F(U(a, £)) is the flux across interface a at 
time t. 


In regions where U is smooth, (4.1) is equivalent 
to the non-conservation form 


U t + A C U X = 0, (4.4) 


where the Jacobian matrix A c is 


A c = 


<9F 

au' 


(4.5) 


With F-C denoting the fc-th component of F, 


where 

«32 = -3(7 - 1)m 2 /2 + 7e//9. 

For the advection equation, the speed is a; here, 
there are three wave speeds; they can be calculated 
by diagonalizing A c . 

Since the centered schemes are very simple, they 
are described first. Also note that we consider only 
the interior points; boundary conditions are stan- 
dard, and are omitted. 

4.1. The N-T and CE/SE schemes. Let 

Xj + i /2 = (j 4- 1/2 )h be the cell interfaces and Xj = 
jh the cell centers, j = 0, 1,2, . . . Recall that the 
reconstruction cell j is the interval [xj^i , Xj + i] of 
length 2 h. 

At time t n , assume that we know Uj (i.e. U”) for 
all j ; Uj appoximates the average of U on the re- 
construction cell j. We wish to calculate U" +1 . In 
the CE/SE case, the slopes {(U x )j} are also stored, 
and {(U x )" +1 } must also be calculated. 

For stability, the time step At is required to sat- 
isfy the CFL condition 

|max(|uj| + aj ) 


At 


< 1, 


(4.7) 


where a denotes the speed of sound: a = ('yp/pY 1 ' 2 - 
Notice that the time step At here corresponds to 
the quantity At/2 in Chang (1995) and Wang and 
Chang (1999). 

To update the cell average U” +1 , the data Uj is 
ignored. The midpoint rule for (4.3) on the recon- 
struction cell [xj-i,Xj+Y\ takes the form 


U »+! = u * + ^( F "+ 1 /2 

3 3 2 h V J- 1 


T7'™+l/2\ 

F j+ 1 )’ 


(4.8) 


where the start-off value U* takes the place of Uj 
and will be defined below. 

Next, we describe the reconstruction function. In 
the CE/SE case, the stored slope is employed; for 
the N-T case, the slope is estimated by Van Al- 
bada’s weighted average: 


( u *)l = X wtav(Uj+i - Uj, Uj - Uj_i). (4.9) 


(A-c)m 


<9F W 
<9U« ' 


By rewriting F in terms of p, m ( m 
after some algebra, 


pu), and e, 


The time derivative is given by (4.4): 

(U t )j = — (A c )j(U x )j. (4.10) 

The reconstruction function takes the form 


A r = 


0 1 0 

(7 — 3)u 2 /2 (3 — 7 )u 7— 1 

(7 — 1 )u 3 — 7«e / p 0,32 7 u 


(4.6) 


=V j + (V x ) j (x-x j )+(U t ) j (t-t n ). (4.11) 

When updating the solution at j , rj is ignored 
while rj_i and rj + i are employed as shown in 
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Fig. 3.1. The start-off value is evaluated at t — t n , 


u* = \ [Uj_! + |(U x )j _ 1 + 
U i+ 1 -*(U*) i+1 ]. 


(4.12) 


To estimate the flux F” +l ' 2 in (4.8) for an arbi- 
trary j, one first obtains 

r j (x j ,t n+1/2 ) = \J j + f(V t ) j , (4.13) 


spatial derivative is estimated by Van Albada’s 
weighted average: 

(V*)i = i wtav(V j+1 - Vj,Vj - Vj—,). (4.20) 
The time derivative then follows from (4.17): 

(Vt)j = —(A P )j (V x )j. (4.21) 

The linear reconstruction takes the form 


and then evaluates the flux at this state. The above 
completes the algorithm of the N-T scheme. 

Instead of (4.13), which was utilized by Nessyahu 
and Tadmor (1990), Chang (1995) employed 

(Ft),- = (Ac)j(U t )j, (4.14a) 

and 

F ; +1/2 = F, + ^(F<), • (4-146) 

This calculation is costlier than (4.13), but it is use- 
ful in the derivation of the reversible scheme. 

Finally, to update the slope for the CE/SE 
scheme, the point value U" +l is calculated by the 
reconstruction, 

U" +1 = Uj + (Ut)jAt, (4.15) 

and 

(U x )? +1 = 4 wtav(U"+ 1 - U” +1 , U” +1 - U"+ x ) . 

(4.16) 

The above completes the description of the CE/SE 
method. 


4.2. Second-order upwind scheme. Here, 
for reasons of economy, we interpolate the primitive 
variable V. Equation (4.4) implies 




V t + A P V X 

= o, 


(4.17) 

where 

( P\ 

( u 

p 

0 \ 


V = 

4 

A p = I 0 

u 

1/p 

. (4.18) 


w 

\o 

IP 

u ) 



Note that A p is simpler than A c . As a result, (4.17) 
is slightly more economical than (4.4). 

By applying the midpoint rule to the integral 
form (4.3) on the cell j, 

Uj +1 = U, + ^(F;* 1 ; 2 - F£$). (4.19) 

The problem reduces to obtaining {Fj^y 2 }. To 
this end, first {V ; } are calculated in a straight- 
forward manner from their definitions. Next, the 


(4.22) 

At each interface j + 1/2 time f” +1 / 2 , as shown in 
Fig. 2.1, the reconstruction rj yields a value de- 
noted by V L , and the reconstruction rj + i , ~V R : 

v l = T j (x j+1/2 ,f n+1/2 ), and 

v r = r i+i(xj+i/2,t n+1/2 )- 

With the above V L and V R , the flux is calculated 
by upwinding as described below. 

4.3. Upwind flux. The upwind flux employed 
here is Roe’s splitting (1986) with an entropy cor- 
rection. Roe’s scheme amounts to picking the up- 
wind side in the characteristic frame depending on 
the wave speed. 

Recall that for the case of advection, the upwind 
flux is 

fu = H au L + au R ) - \\ a \( u R - U L )> 
where a is the wave speed. In other words, with 

Sr~ Jl ~ a ( u R ~ u l)i 

dissipation is added by replacing a with |a|. For the 
Euler equations, A c = <9F/<9U. The question is: 
does the mean value theorem hold for these equa- 
tions? More precisely, with V L and \ R given by 

(4.23), denote AF = F r -Fl and AU = Ur — U^. 
Then, is there a state U which satisfies 

AF = A C AU ? 

If such a state exists, the upwind flux is given by 
Fu= 4(F l +Fr)-||A c |AU, 

where |A C | is calculated via a diagonalization. 

Another way to explain the motivation for Roe’s 
state is as follows. With V L and V R given by 

(4.23), which state should we use for the diagonal- 
ization? We could use the average state \ (Vl + 
V fi ), but when the wave speed is zero, such a state 
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leads to ambiguity in the upwind side selection. 
This ambiguity can be avoided by using Roe’s tilde 
state. 

To diagonalize A c , denote the Jacobian matrix of 
the transformation between the primitive and con- 
servative variables by M (Hirsch 1990): 

Then 

( 1 

M = u 

\u 2 /2 pu 1/(7 - 1 )/ 

By applying the chain rule to (4.4), 

A p = M -1 A C M. (4.26) 

Thus the diagonalization of A c reduces to that of 
A p . The eigenvalues of A v are 

A (1) = u - a, A (2) = «, A (3) = u + a. (4.27) 


M = 


<9U 

av 


0 

p 


(4.24) 


(4.25) 


Note that L c is complicated, but we only need R c 
and L p for the final expression of Roe’s scheme. 
Roe’s tilde state is a state U which satisfies 


AF = A C AU, 

(4.35) 

and 

AU = MAV. 

(4.36) 

Here, (4.35) provides consistency to the picking pro- 
cess, while (4.36) leads to the use of L p , which is 
more economical than L c . Solving the above two 
equations, we obtain 

P = y/pLpR, 

(4.37) 

and with 


Pl = Pl/(pl + p), Pr = 1 - Pl, 

(4.38) 

the other two quantities are given by 


u — Plul + Prur, 

(4.39) 


Let R„ be the matrix of the right eigenvectors of 
A p ; L p , that of the left; then, L p = R" 1 . and 


0 

-p/{2a) 

l/(2a 2 ) ^ 


1 

0 

— 1/a 2 

(4.28) 

0 

p/(2a) 

1/(2 a 2 )) 



( 1 

1 1 ^ 


-p ~ 

= ~ a /P 

0 a/p 

(4.29) 


V « 2 

0 a 2 / 



Denote by A be the diagonal matrix whose diagonal 
entries are A^\ A^ 2 \ and A^ 3 ^. Then 

L P A P R P = A or A p = R P AL P . (4.30a, b) 


H =PlH l +P r H r . (4.40) 

With the above tilde state, we can multiply by 
L c to switch to the characteristic frame: 

Gl = L c Fl, and Gr = L c Fr. (4-41) 

The *-th component of either G l or Gr is chosen 
by the sign of A^ for i = 1,2 and 3. The result is 
an upwind flux in the characteristic frame 

G C / = i(G L + G R )-|sign(A)AG. (4.42) 

Next, we switch back to the regular frame by mul- 
tiplying by R c : F f/ = R C G[/, 


The diagonalization of A c follows from the above 
and (4.26): 

L C A C R C — A or A c — R C AL C , (4.31a, b) 
where 

L c = L P M _1 , R c = MR P . (4.32a, b) 
Let H be the total enthalpy, 

H = (e+p)/p. (4.33) 

Then, (4.32b), (4.29), and (4.25) lead to 

/ ! 1 1 \ 

R c = I u — a u u + a I. (4.34) 

\H — ua u 2 /2 H + ua ) 


F u = !(Fl + F„) - |R C sign(A) L C AF. 

Using (4.35) and (4.31b), one obtains 

F C / = i(F L +F H )-iR c |A|L c AU. (4.43) 

Note that the above expression involves no condi- 
tional statement; i.e., due to (4.35), if A'u = 0 the 
choice of the left or the right value in (4.42) leads 
to the same final upwind flux. Next, due to (4.36), 
L C AU = LpAV. As a result, 

F v = |(Fz, + F r ) - |Rc |A| LpAV. (4.44) 

We now discuss the entropy correction. Set 

AW = LpAV. (4.45) 
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For the first characteristic, the spreading rate can At regions where U is smooth, (5.1) is equivalent 
be estimated by (for details, see (Huynh 1995a)), to the non-conservation form 


AX — AV 1 * — ^ + (4.46a) 

2p 

or if it is the third one, 

AA = AA (3) = (7 + l)«A W (3 )^ 

(4.466) 
2 P 

For any real number z, define the negative part of 
z by z~ = min(^, 0), and the positive parts of z by 
2 + — max(^,0). To see how far the wave speeds 
spread past zero, evaluate 

r, = (— |A| + §AA)+. (4.47) 

For an approximation of Osher’s flux (1984), set 

7 ? =(-|A| + i|AA|)+. (4.48) 

Roe’s splitting with an entropy correction can be 
coded as follows. 

Given the left and right states, calculate Hi and 
H r by (4.33) and the tilde state by (4.37)-(4.40). 
Next, if u > 0, obtain Au/ 1 ) via (4.45), AA^ 
(4.46a), rf 1 ^ (4.47) and, with R). denoting the first 
column of R c , 

Ft; = F l + [(u - a)~ - A^jAw/bRi. (4.49a) 

Otherwise, obtain via (4.45), AA^ (4.46b), 

rj( 3 ) (4.47), and 

Fu=F r - [(fi + a)+ + |t ? ( 3 )]Aw( 3 )r3. (4.496) 


5. Two-dimensional extensions on a 
quadrilateral mesh. The second-order upwind, 
N-T, and CE/SE schemes are extended to the 2D 
case in this and the next section. 

The 2D Euler equations take the form 


U( + A C U* + BpUy — 0; (5-2) 

here, with F ( 0 denoting the fc-th component of F, 


<9F (& ) 

{A c )k,i — (®c)<!,/ 


<9G « 


Let V be the vector of primitive variables. Then 
(5.2) can be put in the form 


V t + A P V X + B p V y = 0. (5.3) 

As in the ID case, A p and B p are simpler than A c 
and B c . 

Let 0 be a spatial control volume whose bound- 
ary is <90. At each point on the boundary, denote 
by n — ( n x ,n y ) the outward unit normal and 

F = (F, G). 


Then the normal flux vector is 


F • n = n x F + n y G. (5.4) 


The Euler equations in integral form for the control 
volume O take the form 


d_ 

dt 


SL 


U dxdy ■ 


F ■ n ds = 0. 


an 


(5.5) 


A discussion on the integral form for a space-time 
domain can be found in, e.g., (Roe 1983). 

Below, O is a polygon of index j with n e edges. 
Denote by area, the area of D, and by U j the aver- 
age value of U on fl at time t n . For each i-th edge, 
denote by its length, fit its outward unit normal, 
and F" +1 / 2 the flux vector F at the midpoint of 
the edge at time f” +1 / 2 , i = 1 ,...,n e . Then, for 
second-order accuracy, (5.5) can be approximated 
by 


where 


U 4 +F(U) x + G(U) y -0, 


U 


/ P 

pu 
pv 
e 


V e / 


(5.1) 


F = 

( P 2 U \ 
I pa + p 

, and G = 

( ^ \ 

pvu 


puv 

V (e + p)uJ 


pv z + p 
V (e + p)v / 


U " +1 (F" +1/2 • *). (5.6) 

Next, let the ( x , |/)-plane be divided into nonover- 
lapping quadrilaterals called cells; two adjacent 
quadrilaterals share a common edge (Fig. 5.1). The 
mesh can be structured: the cells are indexed by i 
and j ; it can be unstructured: the cells are indexed 
by j and, for each cell j, the four neighboring cells 
are located by pointers. Without loss of generality, 
we assume the mesh is unstructured. 
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Fig. 5.1. A quadrilateral cell and its four neighbors. 
The centroids of the cells are denoted by Q, A, B, 
C , and D. For the second-order upwind scheme, to 
update the solution at Q, we need the fluxes at time 
j.n+i /2 midpoints of the cell edges, denoted by 

A', B', C', and D' . 

5.1. Second-order upwind scheme. Suppose 
the data {Uj} (i.e., {U"}) are known; U j approx- 
imates the average of U on the cell j. We wish to 
calculate {U” +1 }. 

Let EFGH be a typical quadrilateral whose cen- 
troid is Q ; this quadrilateral is also identified as the 
cell Q. Let A, B , C, and D be the centroids of the 
four neighboring cells as shown by Fig. 5.1. 

Applying (5.6) to the cell Q, the problem reduces 
to obtaining the fluxes at time t" +1 / 2 across the 
four edges. These fluxes are evaluated at the mid- 
point A 1 , B' . C" , and D' of the edges by using the 
reconstruction functions described below. 

For reasons of economy, we interpolate V. The 
spatial derivatives of V at Q are estimated by Van 
Albada’s weighted average in the following manner. 
Set 

e 4 = i(F^ + Gft) and S n = \(gP + W&). (5.7 a) 
We can also use 

e*£ = \ cA and e v — t,DB. (5.76) 

Next, let e*£ and e v be the basis vectors for the (£, g) 
coordinates: 

(x,y) = £ +g e v . (5.7c) 

Then 

(V € ) Q =wtav(V A -V Q ,V Q -Vc), 

(V,) Q = wtav(Vfl — Vq, Vq — V_d). 

From Vj and V I( , the chain rule yields V x and V v . 
The time derivative follows by applying (5.3): 

(Vt)o = -(A P ) Q (V x ) q - (B p ) q (V„)q. (5.9) 


The linear reconstruction takes the form 
r Q(x,y,t) =V Q + (V x ) Q (x - x Q ) + 

(V v )«,(» -tfo) + (V t )g (*-*")■ 

(5.10) 

Note that the term ‘reconstruction’ is typically used 
for tq(x, y, t n ). Here, we use the term ‘reconstruc- 
tion’ in the extended sense above. 

At time f" +1//2 , at the midpoint of each of the 
four edges, say, at A', the reconstruction tq yields 
a value denoted by V L ; the reconstruction r^, a 
value denoted by V R : 

V L =r Q (x A ' -xq, y A , -y Q , f n+1/2 ), 

V R =r a(x A ’ - x A ,VA’ -VA,t n+1/2 ). 

(5.11) 

With the above V L and V R , the flux is calculated 
by upwinding as described in the next subsection. 
The other three fluxes at B ’ , C . and D’ are similar. 

Observe that for this upwind scheme, there is a 
trade-off between computing time and storage if the 
quadrilateral mesh is unstructured. This observa- 
tion also holds for an unstructured triangular mesh, 
but if the mesh is structured, it does not apply. For 
each cell Q, we can carry out five reconstructions 
(per cell) at Q, A, B. C, and D, and four upwind 
fluxes to update the solution at Q. If the number of 
cells is N , this approach needs to evaluate 5 N (vec- 
tor) reconstructions and 41V upwind fluxes. The 
number of values {U?} and {U” +1 } to be stored 
is 2 x (41V). To reduce computing time, we can 
carry out only one reconstruction per cell and store 
the left and right interface values at time t” +1 / 2 
for all edges, and then, loop over all edges to get 
the fluxes. If the number of cells is IV, the num- 
ber of edges is roughly 21V. This approach needs 
to evaluate only N reconstructions and 2 N upwind 
fluxes. However, the number of additional values to 
be stored is 4 x (41V) . 

5.2. Upwind flux. Again the upwind flux em- 
ployed here is Roe’s splitting (1986) with the en- 
tropy correction in Huynh (1995). 

Let n = ( n x ,n y ) be the outward unit normal 
at A’ . If U is a state at A ', then, because the flux 
is homogeneous of degree one in U, 

F n = (n x ,n y )-( F,G) = (n a ,A c +n y B c )U. (5.12) 

To obtain the upwind flux, we need to diagonal- 
ize the matrix (n x A c + n y B c ). (Notice that be- 
cause the tangential component is ignored, this di- 
agonalization is sometimes said to be dimensional 
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splitting; but the flux for the centered schemes is 
also given by (5.12) where the tangential compo- 
nent is ignored. As such, the upwind and centered 
schemes here are equally dimensional splitting or 
non-splitting.) Next, set 

q = n x u + n y v. (5.13) 

Then the eigenvalues of the matrix (n x A c -\-nyB c ) 
are 

A ^ = q — a, A ^ = q, A^ = q, \^=q + a; 

(5.14) 

and, the matrix R c of the right eigenvectors is 

/ 1 0 1 1 \ 

u — n x a —n y a u u + n x a 

v — n y a n x a v v + n y a 

\ H — qa a(-n y u + n x v ) “ \ v H + qa ) 

(5.15) 

As for the matrix of the left eigenvectors of n x A p + 
n y TH p , where the subsript p stands for ‘primitive’, 


f° 

-n x p/(2a) 

-n y p/(2a) 

l/(2a 2 ) 

0 

~n y p/a 

n x p/a 

0 

1 

0 

0 

-1/a 2 

\0 

n x p/(2a) 

n y p/{2a) 

1/(2 a 2 ) 


(5.16) 

We only need the first and last columns of (5.15) 
and the first and last rows of (5.16). 

Roe’s splitting with an entropy correction can be 
coded as follows. 

Let the left and right states and the unit nor- 
mal n be given. First, calculate Hl and Hu by an 
expression similar to (4.33) and the tilde state by 
(4.37)-(4.40) (v is similar to u and H). Next, if 
q > 0, obtain An/ 1 ) via (4.45), AA^ (4.46a), t/ 1 ) 
(4.47) and, with R). denoting the first column of 
R-c 

(F • n)u = (F • n) L + [(q - a)~ - Au> ( 1 ) Rj. 

(5.17a) 

Otherwise, obtain Aw^ via (4.45), AA^ (4.46b), 
q (4.47), and 

(F • n) v = (F • n) R - [(q + d)+ + |») ( 4 )]Aw( 4 )R<. 

(5.176) 

5.3. Second-order centered schemes. The 

CE/SE schemes were extended to the case of a 
quadrilateral mesh by Zhang et al. (2002). This 
approach to extension is also applied to the N-T 
scheme below. The key difference here, however, 
is that the slope estimates are simplified, and they 
only use the data at the current time level. The 



Fig. 5.2. Reconstruction cells for the two centered 
schemes. The domain where the reconstruction is 
valid at time t = t n for the cell Q is the dotted line 
octagon, namely, AEBFCGDH . This octagon is 
called the reconstruction cell Q here. It is also the 
control volume on which fluxes are balanced when 
the solution at Q is updated. The centroid of the oc- 
tagon is marked by the gray dot and denoted by Q . 

simplified slope estimate together with an observa- 
tion on distributing the fluxes to the neighboring 
cells result in an extended N-T scheme which is 
faster and requires considerably less storage than 
the CE/SE method for the case of an unstructured 
mesh. The extended CE/SE scheme here is also 
different from that by Zhang et al. (2002) in that 
essentially, the scheme which carries along the in- 
terface values explained after (3.14) is employed. 
Finally, the presentation below is simpler. 

Consider a typical quadrilateral EFGH whose 
centroid is Q shown in Fig. 5.2. Let the centroid of 
the four neighboring cells be denoted by A , B , C, 
and D. For the two centered schemes, let the re- 
construction cell Q be defined as the (dotted line) 
octagon AEBFCGDH . At time t = t rl , the re- 
construction for these two schemes is assumed to 
be valid on this octagon. This octagon is roughly 
twice as big as the original cell EFGH. The cen- 
troid of this octagon is marked by a gray dot and 
denoted by Q. It is called a solution point by Wang 
and Chang (1999). Note that the two adjacent re- 
construction cells Q and A overlap on the quadri- 
lateral QHAE (see also Fig. 5.3). Also note that 
the reconstruction cell is the spatial part (or spatial 
projection) of the solution element in Zhang et al. 
( 2002 ). 

At time t n , assume that we know Uj (i.e. U") 
for all j; Uj appoximates the average of U on the 
reconstruction cell j. We wish to calculate UJ 1 1 . 

If Uq approximates the average of U on the re- 
construction cell Q, it is considered to be the value 
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U at the centroid Q, not at Q. A more precise no- 
tation would be Uq, but for simplicity, we use the 
notation Uq below. 

In the CE/SE case, the spatial slopes {(U^} 
and {(Uj,)^} are also stored; in addition, {(U x )" +1 } 
and {(U y )" +1 } must be calculated. 

To update the cell average Uq + 1 . the data Uq 
is ignored. Denote by area(rc(Q)) the area of the 
reconstruction cell Q. Then the midpoint rule for 
(5.5) on this reconstruction cell takes the form 

area(rc(Q)) Uq + 1 = area(rc(Q)) Uq — 

Ai E- = 1 ^(F ” +1/2 •*), 

(5.18) 

where the start-off value Uq and the eight fluxes 
are defined below. 

First, we describe the reconstruction function for 
an arbitrary cell, say, cell Q shown in Fig. 5.2. In 
the CE/SE case, the stored slopes are employed; 
for the N-T case, they are estimated by using the 
neighboring values as follows. Let the (£,i?) coor- 
dinates be defined as in (5.7c) with = \CA and 
e v — \DB, or we can also use (5.7a). Then 

(U ? )q = wtav(U+ - Uq,Uq - U c ), 

(5.19) 

(U„)q = wtav(U B - Uq,Uq - U d ). 

From U^ and U tj . the chain rule yields U J; and U y . 
The time derivative follows by applying (5.2): 

(Uj)q = -(A c )q (U x )q - (B c )q (U,)q. (5.20) 

Since Uq approximates the average of U on the 
reconstruction cell, it is considered to be the value 
at the centroid Q. The linear reconstruction takes 
the form 

T Q {x,y,t) = Uq + (U^q (x -xq) + 

(U v )q ( y - yq) + (U t ) Q {t - t n ). 

( 5 - 21 ) 

Note xq and yq above take the place of xq and yQ 
in the upwind case. 

Next, when updating the solution at Q , tq is 
ignored, and r a, r/j. re, and r o are employed 
(Fig. 5.3). The start-off (average) value at time 
t = t n is calculated as follows. Denote by A* the 
centroid of QHAE ; B * , QEBF; C * , QFCG; and, 
D *, QGDH. At time t = t n , the value at A* is 
evaluated using r^; at B * , using r B ; C*, tq; and 



Fig. 5.3. Centered schemes. When updating the so- 
lution at Q, vq is ignored, while the reconstructions 
at the neighboring cells r r B , vq, and are em- 
ployed. Denote the centroid of QHAE by A*, and 
that of QEBF, B* . The values needed from are: 
at time t n at A*; at time t" +1 7 2 at the midpoint of 
AH and the midpoint of AE marked by (+)■ The 
values needed from r B are: at time t n atB*; at time 
£n+i /2 a £ £^ e misprint 0 f BE and the midpoint of 
BF marked by (+). Similar statements hold for the 
reconstructions at C and D. 

D*, r B - The start-off value Uq is given by 

area(rc(Q)) Uq = 

area(QiLAE) r+(:rA. ,yA*,t n ) + 
area(Q£’5F) r B (x B . ,i/ B .,t”) + 
area(QFCG) rc{xc* , yc * , f") + 
area (QGDH) r D (x D . , y D * , t n )\ 

(5.22) 

here, again, area(rc(Q)) is the area of the recon- 
struction cell (octagon) AEBFCGDH . 

As for the fluxes j?” +1 / 2 . across the eight 
edges, the fluxes across AH and AE are calculated 
by the reconstruction r ,i ; those across BE and BF, 
by r B ; CF and CG , by tq; and DG and DH, by 
td- Denote by M the midpoint of AE. Then one 
calculates the conservative variables at M at time 
^n+i /2 v - a r a{%m ,VM, t n+1 / 2 ), and the flux across 
AE by (5.4). The other seven fluxes are obtained in 
the same manner. This completes the quadrilateral- 
mesh extension of the N-T scheme. 

Next, the slope update for the CE/SE scheme be- 
low is a simplification of that in Zhang et al. (2002). 
First, the point value at A at time t n+1 , denoted by 
U l +1 , is calculated by the reconstruction, 

U(t +1 = r A (x A ,y A ,t n+1 ) . (5.23) 

The other point values U^ +1 , Uq + 1 , and U^ +1 are 
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H 


D 

Fig. 5.4. Distributing the fluxes results in a 
faster algorithm with no penalty in storage for 
the extended N-T scheme. After calculating 
V x , U y , and U( at Q, distribute the quantity 
axea,(QHAE)rQ(xA*,yA'‘,t n ) and the two fluxes 
at H' and EJ to cell A; distribute the quantity 
axea,(QEBF) tq{xb* , t/s* , t n ) and the two fluxes at 
E' and F' to cell B; the distribution to cells C and 
D is similar. 

obtained in a similar manner, and 

(U 4 )q = wtav(Uf 1 - Uq +1 , Uq +1 - U£ +1 ), 

(Ujj)q = wtav(U£ +1 - Uq +1 , Uq +1 - U£ +1 ). 

(5-24) 

The spatial derivatives (Ua,)g +1 and (U^)^ 1 can 
then be calculated via the chain rule. This calcula- 
tion completes the quadrilateral-mesh extension of 
the CE/SE method. 

The following observations are in order. 

(a) Contrary to the upwind case, for an unstruc- 
tured quadrilateral (or triangular) mesh, the above 
extended N-T scheme can be coded in a manner 
which results in a faster algorithm with no penalty 
in storage. Indeed, instead of gathering fluxes to 
update the solution, fluxes can be distributed to 
the neighboring cells in the following manner. 

First, set all U" +l to zero. Next, evaluate 
area(rc(j)) U” 1 1 for all j as follows. In a Do loop 
with index j, when j equals Q , obtain tq. Calculate 
the values of tq at time t = t n at A*, B*, C*, and 
D* (see Fig. 5.4). Calculate the values of tq at time 
t = t n+ !/ 2 at the spatial locations marked by (+) 
in Fig. 5.4. These locations are: E ' , the midpoint 
of QE\ Ff that of QF\ G', QG\ and Hf QH. Dis- 
tribute the quantity area(QJJ'A.E) , t/A* ,f") 

and the two fluxes at H' and E' to cell A. and 
store the sum in U(^ +1 ; distribute the quantity 
area (QEBF) ta(xb* , 2 /b*, f") and the two fluxes at 
E 1 and F' to cell B , and store the sum in U^ +1 ; the 
distribution to cells G and D is similar. When the 



Fig. 5.5. Two staggered meshes in Arminjon et 
al. (1995) overlay to form a nonstaggered mesh in 
Zhang et al. (2002) provided that the indices ( k,l ) 
are assigned along the two diagonal directions. 

loop is completed, each cell will have received all 
the fluxes it needs to update its cell average value. 

(b) The above extension of the nonstaggered N-T 
scheme relates to the extension of the staggered ver- 
sion in Arminjon et al. (1995) and Jiang and Tad- 
mor (1998) as follows. 

For the staggered version, the scheme alternates 
between the black and gray dots shown in Fig. 5.5. 
More precisely, let a mesh of squares whose edges 
are the black lines, centers black dots, and vertices 
gray dots be given. At time t n , the values U"^ 
are known at the black dots. The reconstruction 
at each black dot is assumed to be valid on the 
corresponding square (black lines). At time t n+1 , 
calculate the solutions at the gray dots by balancing 
fluxes on the squares whose edges are the gray lines. 
At time t n+ ' 2 , obtain the solutions at the black dots 
by balancing fluxes on the squares whose edges are 
black lines. 

This extension of the staggered scheme can be 
made nonstaggered in the following manner. At 
time t n , assume the data are known at both the 
black and gray dots. At time t n+1 , obtain the 
solutions at black dots by using the reconstruc- 
tions at gray dots, and the solutions at gray dots 
by using the reconstructions at black dots. If we 
assign indices (k,l) along the two diagonal direc- 
tions (observed by Drs. Ananda Himansu and Xiao- 
Yen Wang), the resulting scheme is very similar to 
the nonstaggered extension of the N-T scheme pre- 
sented here. 

(c) Compared with the CE/SE scheme discussed 
in Zhang et al. (2002), we can save storage by 
carrying along the point values Ug instead of the 
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Fig. 6.1. The centroid of the triangle EFG is de- 
noted by Q. The centroids of the three neighboring 
cells are denoted by A, B, and C. For the second- 
order upwind scheme, to update the solution at Q, 
we need the fluxes at time t n+ x ' 2 at the midpoints 
of the cell edges, denoted by A' , B' , and C' . 


i 



Fig. 6.2. To derive the triangular-mesh version for 
the weighted average formula, it is convenient to 
denote A, B, C by 1, 2, and 3 respectively. 

slopes. In addition, we can save computing time by 
applying the observation on distributing the fluxes 
so that each cell is visited only once instead of four 
times. Also note that for a quadrilateral mesh, since 
there is no reversible scheme, it is not clear how to 
adjust numerical dissipation even if the slopes are 
carried along. 

(d) For the CE/SE method, the flux at time 
t n+ 1/ 2 can be obtained via (Fj)^ and (G t )A in a 
manner similar to (4.14). Such a calculation, how- 
ever, is costlier. 

6. Two-dimensional extensions on a trian- 
gular mesh. For the schemes discussed here, the 
extension to a triangular mesh differs from that to a 
quadrilateral mesh described in the previous section 
only in the slope evaluation, which employs three 
neighboring data, not four. 

Let the (x,y)-p\ane be divided into nonoverlap- 
ping triangles called cells; two adjacent triangles 


share a common edge (Fig. 6.1). Assume that the 
mesh is unstructured: the cells are indexed by j 
and, for each cell j, the three neighboring cells are 
located by pointers. 

Let EFG be a typical triangle whose centroid is 
denoted by Q. This triangle is also identified as the 
cell Q. Let A, B, and C be the centroids of the 
three neighboring cells as shown in Fig. 6.1. 

Next, we describe the triangular-mesh extension 
of Van Albada’s weighted average. The following 
version is a modification of the extension by Wang 
and Chang (1999). It is designed for the upwind 
and the N-T schemes where only data at the cur- 
rent time level are available. It also makes the ob- 
servation on distributing the fluxes applicable. 

For the purpose of deriving this weighted average, 
as will be seen in (6.1), it is more convenient to use 
numbers 1, 2, and 3 to denote the neighbors; e.g., 
instead of (xa,Va), we use (*1,2/1) (see Fig. 6.2). 
Let uq be a scalar value at the the point (xq,j/q), 
and similarly, u%, M2 , and «3 , at the three neighbor- 
ing cell centers. Next, the values at any three points 
determine a plane in the ( x,y,u ) space; for exam- 
ple, the values Mi, U2, and M3 at (a?i , 3/1), (^2^2)5 
and (*3,1/3) respectively, determine a plane. Denote 
the slopes of this plane by (m x )i 23 and (m j/ )i 23 - 

The plane we wish to obtain is biased toward the 
least steep one among the three planes 12 Q, 23 Q, 
and 31 Q, but if ( xq,pq ) lies on one of the edges 
of the triangle 123, then one of the three planes is 
not well defined. To avoid this problem, let (*o,2/o) 
be the centroid of the triangle 123. The value Mo is 
defined by linear extrapolation using the plane 123: 

Uq = UQ + (u x )i23(X0 ~ %q) + (Wj/)l23(l/0 ~ Vq)- 

The final slopes are obtained by the values at the 
points 0, 1, 2, and 3 as follows. After calculating 
the slopes of the three planes 012, 023, and 031, set 

#1 = [(Wz)o23] 2 + [(%)o23] 2 i 

02 — [(m x )o 3 i ] 2 + [(%)o3l] 2 , (6.1) 

03 = [(«:c)oi 2 ] 2 + [(M y )oi 2 ] 2 , 

and 

0123 = 0102 + 0203 + 0301- 

Note the less steep the plane 012, the smaller the 
quantity 03 . The final slope at Q is given by biasing 
toward the least steep plane among the above three: 

(u x )q = [(0102)(«s)oi2 + (0203)(«;c) 023 + 
(030l)(Mi)o3l]/0123- 

(6.2) 
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F 


Fig. 6.3. Reconstruction cells for the two centered 
schemes. The domain where the reconstruction is 
valid at time t = t n for the cell Q is the dotted 
line hexagon , namely , AEBFCG. This hexagon is 
called the reconstruction cellQ. It is also the control 
volume on which fluxes are balanced when the solu- 
tion at Q is updated. The centroid of the hexagon 
is marked by the gray dot and denoted by Q. 

A similar expression holds for ( u v )q . This 

completes the triangular-mesh extension of the 
weighted average. 

6.1. Second-order upwind scheme. Suppose 
the data {Uj} (i.e., {U"}) are known; U ; - approx- 
imates the average of U on the cell j. We wish to 
calculate {U" +1 }. 

To update the solution at the cell Q, we need to 
calculate the fluxes at the midpoint A ', B', and C' 
at time t n+1 / 2 (Fig. 6.1). These fluxes are obtained 
by using the reconstructions described below. 

Consider an arbitrary cell, say, cell Q. The values 
V at Q and at the three neighboring cells determine 
the slopes and V y at Q via (6.2). The time 
derivative then follows from (5.9), and the recon- 
struction, (5.10). At time t n+ 1 /2 , at say, A\ the re- 
construction vq yields V, . and rx yields V R . The 
flux at A! is calculated by the upwind step (5.17). 
The other two fluxes are similar. These fluxes com- 
plete the upwind algorithm. 

Observe that for this upwind scheme, there is a 
trade-off between computing time and storage. For 
each cell Q, we can carry out four reconstructions 
(per cell) at Q, A, B , and C, and three upwind 
fluxes to update the solution at Q. To reduce com- 
puting time, we can carry out only one reconstruc- 
tion per cell and store the left and right interface 
values at time t" +1 / 2 for all edges, and then, loop 
over all edges to get the fluxes. 

6.2. Second-order centered schemes. The 

CE/SE schemes were extended to the case of an 


unstructured triangular mesh by Wang and Chang 
(1999). This approach to extension is also applied 
to the N-T scheme below. The key difference here, 
however, is that the slope estimates are modified 
so that only the data at the current time level 
are needed. As in the quadrilateral-mesh case, the 
modified slope estimate made the observation on 
distributing the fluxes to the neighboring cells ap- 
plicable. This observation, in turn, results in an 
extended N-T scheme which is faster and requires 
considerably less storage than the CE/SE method 
(slopes need not be stored). The extended CE/SE 
scheme here is also different from that by Wang 
and Chang (1999) in that essentially, the scheme 
which carries along the interface values explained 
after (3.14) is employed. Finally, the presentation 
below is simpler. 

For the two centered schemes, let the reconstruc- 
tion cell Q be defined as the (dotted line) hexagon 
AEBFCG (Fig. 6.3). At time t = f", the recon- 
struction for these two schemes is assumed to be 
valid on this reconstruction cell. The centroid of 
the reconstruction cell is marked by a gray dot and 
denoted by Q. The two adjacent reconstruction 
cells Q and A overlap on the quadrilateral QGAE 
(Fig. 6.4). Note that the reconstruction cell is the 
spatial part (or spatial projection) of the solution 
element in Wang and Chang (1999). 

At time t n , assume that we know Uj (i.e. U") 
for all j ; Uj appoximates the average of U on the 
reconstruction cell j. We wish to calculate U” +1 . 

In the CE/SE case, the spatial slopes {(tXOj) 
and {(U y )j} are also stored; in addition, {(U x )" +1 ) 
and {(U y )” +1 } must also be calculated. 

Next, we describe the reconstruction function for 
an arbitrary cell, say, cell Q shown in Fig. 6.3. In 
the CE/SE case, the stored spatial slopes are em- 
ployed; for the N-T case, the value Uq, which is 
considered to be at Q , and those at A, B. and C 
determined U T and U y via the weighted average 
(6.2). The time derivative is given by (5.20), and 
the reconstruction, (5.21). 

The solution at Q is calculated by balancing 
fluxes over the reconstruction cell Q. Denote by 
area(rc(Q)) the area of the reconstruction cell Q. 
Then, similar to (5.18), by the midpoint rule, 

area(rc(Q)) Uq + 1 = area(rc(Q)) Uq — 

At EL h(.F ? +1/2 • ft *). 

(6.3) 

When updating the solution at Q, tq is ignored, 
and rx, r#, and tq are employed (Fig. 6.4). Denote 
by A* the centroid of QGAE\ B * , QEBF; and C* , 
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Fig. 6.5. For the staggered version, the scheme can 
alternate between any two of the three sets of data: 
at cell vertices , at cell centers , and at midpoints of 
edges. 


Fig. 6.4. Centered schemes. When updating the 
solution at Q, tq is ignored, while the reconstruc- 
tions at the neighboring cells r A , rg, and r c are 
employed. Denote the centroid of QGAE by A*; 
that of QEBF, B* ; and that of QFCG, C* . The 
values needed from r a are: at time t n at A*; at time 
j.n+i /2 midpoint 0 f AE an d the midpoint of 

AG marked by (+). The values needed from r B are: 
at time t n at B* ; at time t n+1 ^ at the midpoint of 
BE and the midpoint of BF marked by (+). A sim- 
ilar statement holds for the reconstruction at C . 

QFCG. The start-off (average) value is given by 
&re&(AEBFCG) Ug = 

axea,(QGAE)r A (x A »,y A ^,t n ) + 

(6.4) 

a,ma,(QEBF)r B (x B *,yB',t n ) + 

&re&(QFCG)r c (xc*,yc*,t n ) ■ 

As for the fluxes F/ +1// ' 2 . n i across the six edges 
at time t" +1 / 2 , those across AG and AE are cal- 
culated by the reconstruction r^; those across BE 
and BF, by r^; and CF and CG, by r c- Denote 
by M the midpoint of AE. Then one calculates the 
conservative variables at M viar j 4 (xM,t/M,f ra+1 / 2 ), 
and the flux across AE by (5.4). The other five 
fluxes are calculated in the same manner. This 
completes the triangular-mesh extension of the N-T 
scheme. 

Next, the slope update for the CE/SE scheme be- 
low is a simplification of that by Wang and Chang 
(1999). First, the point value at A at time t n+1 , 
denoted by U'( 1 1 . is calculated via the reconstruc- 
tion, 

V n A +1 =r A (x A ,y A ,t n+1 ) . (6.5) 

The point values U^ +1 and Ug +1 are obtained in 
a similar manner. These three values and U g +1 


(at Q) determined (U x )g +1 and (UjJg -1-1 via the 
weighted average (6.2). This calculation com- 
pletes the triangular-mesh extension of the CE/SE 
method. 

The following remarks are in order. 

(a) Distributing the fluxes results in a faster al- 
gorithm with no penalty in storage for the ex- 
tended N-T scheme. Indeed, let E' denote the 
midpoint of QE; F', QF\ and G', QG. After 
calculating ILj, U y , and U ( at Q. distribute the 
quantity area (QGAE) tq ( x A * , y A * , t n ) and the two 
fluxes at G' and E' to cell A; distribute the quantity 
area (QEBF) tq(xb * , t/s* , t n ) and the two fluxes at 
E 1 and F' to cell B\ finally, distribute the quantity 
axe&(QFCG) tq(xc* , yc* > t n ) and the two fluxes at 
F' and G' to cell C. 

(b) The above observation can also be applied to 
the CE/SE scheme: instead of storing the slopes, 
we can store the point values Ug and, when visiting 
the cell Q, after distributing the fluxes, we update 

(c) The next remark concerns the extension of 
the staggered version. Let an unstructured trian- 
gular mesh of N vertices be given. Then the num- 
ber of cells is roughly 2 N, and the number of edges, 
roughly 3 N. The data can be stored at the vertices, 
the cell centers, or the midpoint of the edges. We 
can alternate between any two of these three sets 
of data, say, between data at vertices and those at 
edges. Then, at one time level, the flow field is re- 
solved by N points (vertices), and at the next time 
level, by 3 N points (edges). Such a solution can 
only be as accurate as resolving the flow field by 
N points only; thus, the scheme is not optimal, 
and there is no gain in accuracy compared with 
the nonstaggered version (see below). In addition, 
the reconstruction, the solution procedure, and the 
bookkeeping for such a method are quite involved. 
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The triangular-mesh extension of the N-T scheme 
here has 2 N pieces of data at each time level, but 
because the reconstruction cell is twice as large as 
the original triangular cell, the flow field is effec- 
tively resolved by N points (not 2 N). However, the 
scheme is simple and has numerous nice features as 
will be shown in the next section. 

7. Fourier analysis for the 2D case. For 

convenience, the 2D extension of the N-T scheme 
presented above will be called the centered scheme 
from here on — as opposed to the upwind scheme. 
It is indeed the centered counterpart of the upwind 
scheme. When more details are needed, we refer 
to it as the Lax-Friedrichs type second-order accu- 
rate centered scheme. It can also be considered as a 
coupled version of the CE/SE method. (Since the 
upwind scheme is sometimes called MUSCL-Roe, 
this scheme could have been named LF-MUSCL- 
NT-CE/SE, but such a name seems too unwieldy.) 
If the slopes are carried along (e = 1/2), the corre- 
sponding 2D extension is called the CE/SE scheme 
with the contributions by MTJSCL, N-T, and the 
various authors explained in the previous sections 
understood. Thus, in this section, we carry out the 
stability and accuracy analyses for the upwind, cen- 
tered, and CE/SE schemes. 

For Fourier analysis, consider the 2D scalar ad- 
vection equation: 

u t + au x + bu v = 0. (7.1) 

The mesh consists of square cells with cell widths 
Ax = Ay = 1 and cell centers at (x, y) = (i,j). 
Next, since Ax = Ay = 1, set 

a x = aAt, and a y = bAt. (7.2) 

For each time step, the data are advected by a spa- 
tial vector (cr x ,a y ), which is called the displacement 
vector (per time step) here. 

In this section, as is typical for Fourier analysis, 
the slope employed is the average slope. 

7.1 Fourier analysis for a square mesh. As- 
sume that a > 0 and b > 0. Then the upwind 
scheme reproduces the exact solution when (a x ^a y ) 
equals (0,0), (1,0), and (0, 1). Loosely put, as can 
be seen from Fig. 7.1, when (< 7 x ,a y ) = (1,0) the 
square centered at B slides onto the square cen- 
tered at Q. 

The centered and the CE/SE schemes, on the 
other hand, do not reproduce the exact solution for 
any (cr x ,cr y ) even when ( cr x ,(j y ) = (0,0). Indeed, 
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Fig. 7.1. The upwind scheme reproduces the exact 
solution when (a x ,(Ty) equals (0,0), (1,0), or (0,1). 
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Fig. 7.2. The centered and the CE/SE schemes do 
not reproduce the exact solution for any ( cr x ,a v ). 

suppose (<j x ,<j y ) — — (1,0) (Fig. 7.2). Then, 

when updating the solution for the cell Q , at the 
beginning of the time step, the contribution from 
the cell B is the reconstruction on BKQE shown 
in Fig. 7.2. By the end of the time step, the re- 
construction on BEGHIK flows into the cell Q via 
the fluxes across BE and BK. Thus, when up- 
dating the solution at Q, the contribution from the 
cell B is the reconstruction on QEGHIKQ, which 
is not the whole reconstruction cell. Therefore, we 
do not recover the value at B. Another way to see 
this fact is by setting all slopes to be zero. Using 
indices i and j and assuming (a x ,a y ) = (1,0), one 
obtains, after some algebra, 

Woj ) 1 = j( 3 u_i,o + uo,-i + uo,i — 1*1,0)- 

That is, we do not recover u- i,o- 

The fact that the solution recovers the exact 
value, loosely speaking, helps keep the error from 
becoming too big. For a triangular mesh, the situa- 
tion is reversed as will be shown in the next subsec- 
tion: the centered and CE/SE schemes recover the 
exact solution for certain values of (a x ,a y ) while 
the upwind scheme does not. 

The upwind scheme yields the following solution 
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at i — 0 and j — 0: 

u ofi l = u o,o + 

\ &x{ — M-2,0 + 5m_i : o — 3uo,o — M l,o) + 

3 °v( — u o ,-2 + 5ito,_i — 3 mo,o — M o,i) + 

3 <7 k( w -2,0 — u -l,0 ~ «0,0 + «l,o) + 

3 &y(UO,-2 — Mo,-l — Uo,0 + Mop) + 

3 (2m_i,_i - m_i, 0 - U - 1,1 - 

Mo,-l + Mo,l — Ml _1 + Ml,o)- 

(7.3) 

The amplification factor A is obtained by replac- 
ing mj in the above expression by e 1 (’ iw *+i w y) , 
where w x and w y are the wave numbers in the x 
and y directions respectively, and 7 = \/—l. The 
exact amplification factor is e 1 (~ <T * w *- <T v w v ) . 

The derivation of the amplification factors for the 
CE/SE scheme is similar to that of the ID case 
except for the more involved algebra. If we carry 
along the interface values, we need to calculate the 
eigenvalues of a 2 x 2 matrix; if we carry along the 
slopes, we need to calculate the eigenvalues of a 3x3 
matrix, but one of the three eigenvalues turns out 
to be identically zero. 

The expressions for the amplification factors of 
the three schemes are lengthy and are omitted. 
However, the stability regions and the plots on ac- 
curacy comparison, which are generated by Mathe- 
matica, will be shown. 

Figure 7.3 shows the stability regions for the up- 
wind, centered, and CE/SE schemes. Note that the 
upwind scheme has the largest stability region, and 
the CE/SE, the smallest. Here, the two amplifica- 
tion factors of the CE/SE scheme yield the same 
stability region. 

To compare errors among the three schemes, first, 
observe that the amplification factors are functions 
of a x , <7 y . w x . and w y . We wish to plot the errors 
for relatively smooth data; therefore, we fix small 
wave numbers w x and w v and plot the errors as 
functions of 

<7 = \(<Tx,<7y)\ 

along the rays a y = constant * o x . (That is, for a 
fixed angle a, a x = acos(a) and a y — asin(a).) 

The phase error per time step is Arg(.4) + (<T x w; x -l- 
cr y w y ). For the second-order schemes discussed 
here, this quantity is proportional to 0(w x )+0(Wy) 
and is the leading error. The dissipation error per 
time step is |.4| — 1, which is proportional to 0(w x )+ 
0(Wy). The total error is \A — e 1 ^ rTxWx ^ rTyWy ^\. 



Fig. 7.3. Stability regions for a rectangular mesh. If 
the displacement vector based at the origin lies in- 
side the stability region, the corresponding scheme 
is stable. The upwind scheme has the largest stabil- 
ity region, and the CE/SE, the smallest. 

For the next four plots, the wave numbers are 
w x = 7t/ 16 and w y = 7t/ 8. The results by the 
standard one-step Lax-Wendroff scheme are also in- 
cluded for ease of reference. Note that this scheme 
has difficulties dealing with shocks. 

In Figs. 7.4, the flow is along the diagonal y = x, 
i.e., o x = (t y . Fig. 7.4(a) shows the phase errors 
as functions of a (denoted by CFL on the plot). 
Note that the phase errors of the centered and the 
CE/SE schemes are nearly identical. 

Figure 7.4(b) shows the dissipation errors as func- 
tions of a. The upwind scheme has the smallest dis- 
sipation error, and the centered scheme, the largest. 

Figure 7.4(c) shows the total errors as functions 
of a. Here, the centered and CE/SE schemes have 
essentially the same errors while the upwind scheme 
has the smallest error. 

Figure 7.5 shows the total errors averaged over all 
flow directions; again, w x — 7t/ 16 and w y = n/8. 

The next four plots in Figs. 7.6 and 7.7 are iden- 
tical to the previous four except that w x = 7t/ 32 
and w y = 7t/16. 

Figures 7.5 and 7.7 show that the total error of 
upwind scheme is considerably less than those of 
the centered and CE/SE schemes, and the latter 
two have essentially identical total errors. 

Note that the order of accuracy of the schemes 
can be verified using these figures. For example, 
the maximum phase error for the CE/SE scheme 
in Fig. 7.4(a) is about 0.0016; that in Fig. 7.6(a) 
is about 0.0002, which is 1/8 of 0.0016. In other 
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Phase Errors; flow along y=x 



Phase Errors; flow along y=x 



CFL 


Dissipation Errors; flow along y=x 



CFL 


Dissipation Errors; flow along y=x 


Errors; flow along y=x 



Fig. 7.4. Errors as functions of a (denoted by CFL); 
flow along y = x; w x = n/16 and w y = tt/8. 

Errors; average all directions 



CFL 


Fig, 7.5. Total errors; average for all flow direc- 
tions; again w x = tt/ 16 and w y = n/8. 



Errors; flow along y=x 



Fig. 7.6. Errors as functions of a; flow along y = x 
w x — 7r/32 and w y = 7r/16. 

Errors; average all directions 



Fig. 7.7. Total errors; average for all flow direc 
tions; again w x = 7r/32 and w y = 7r/16. 
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Fig. 7.8. The upwind scheme does not reproduce 
the exact solution for any (cr x ,o y ) other than the 
obvious value o/(0,0). 



Fig. 7.10. The numbering of the triangles for 
Fourier analysis. 



Fig. 7.9. The centered and the CE/SE schemes re- 
produce the exact solution if (<r x ,o y ) equals one of 
the three vectors A$, B$, or C<$. 


words, when we double the number of mesh points 
in each direction (quadruple the total number of 
points), the wave numbers reduce by a factor of 2, 
the phase error after one time step reduces by a 
factor of 8, and the dissipation error, a factor of 16 
(as can be seen in Figs. 7.4(b) and 7.6(b)). Similar 
observations also hold for the other schemes. 

Note that the errors for the square-mesh case here 
are similar to those for the ID case in §3. 

7.2 Fourier analysis on a triangular mesh. 

Here, each square in Fig. 7.1 is cut into two trian- 
gles along the diagonal from the northwest to the 
southeast corners as shown in Fig. 7.8. 

The first question is when do the schemes recover 
the exact solution? Here, the upwind scheme re- 
covers the exact solution only for the obvious case 
of (a x ,<T y ) = (0,0). The reason is that under the 
translation A($, the triangle A does not match with 
the triangle Q as can be seen from Fig. 7.8. Similar 


observations hold for the translations B('} and C<$. 

The centered and CE/SE schemes, on the other 
hand, do not recover the exact solution when 
(cTx,&y) = (0,0), but they recover the exact so- 
lution when (a x ,<j y ) equals one of the three dis- 
placement vectors A$, FUf), or ct}, as shown in 
Fig. 7.9. Indeed, suppose ( a x ,a y ) = A$. Then, 
loosely put, the reconstruction cell A slides onto 
the reconstruction cell Q. More precisely, at the 
end of the time step, through the start-off value on 
AEQG and the two fluxes across AE and AG, the 
contribution from the cell A is the reconstruction 
on the whole reconstruction cell A. Therefore, af- 
ter one time step, the solution at Q is identical to 
the data at A: we recover the exact solution. This 
property of recovering the exact solution helps the 
centered and the CE/SE schemes gain back some 
accuracy compared with the upwind scheme as will 
be shown later. 

The next few examples using the first-order 
schemes are simple, but they convey the behavior 
of schemes on a triangular mesh. 

For first-order accuracy, the slopes are set to zero, 
and both the centered and CE/SE schemes reduce 
to the L-F scheme. If the data are on a plane, this 
scheme recovers the exact solution. Indeed, the so- 
lution at Q is obtained by applying the displace- 
ment to the plane through the data at A, B and C 
shown in Fig. 7.9. 

The upwind scheme, on the other hand, does not 
recover the plane. Indeed, suppose the flow angle 
is between 0 and — 7t/4 so that the upwind cell to 
cell Q is cell B (Fig. 7.8). Then, when updating the 
solution at Q , the solution involves the data at only 
Q and B. These two pieces of data cannot recover 
a plane (we need three). What happens is that 
typically there is a cancellation of errors (for the 
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fluxes) so that a scheme using a piecewise constant 
reconstruction reproduces the exact solution if the 
data are on a plane. Here, the cancellation of errors 
did not occur. 

The next question is: what does the piecewise- 
constant upwind scheme do to a data which is on 
a plane? It turns out that the scheme propagates 
the plane and turns it into a ‘bumpy plane’. As 
an example, suppose a y = 0 and a x > 0. Then 
(see Fig. 7.10), for the piecewise-constant upwind 
scheme, 

U L,i,j = u L,i,j + 2&x( u R,i—l,j ~ u L,i,j ), (7.4a) 

and 

C — 'RR^iJ 4“ 2a x (uL,i,j R'R^ij)' ( ( .46) 

If the data are on the plane u = x, then (see 
Fig. 7.10), 

u L ,i,j = * - 1 /6, and u RtiJ = i + 1/6. (7.5) 

Next, consider the following data of odd-even type 
noise 

e L,i,j = -1/12) and e Ri j = 1/12. (7.6) 

Suppose the data is a ‘bumpy plane’ obtained by 
superimposing (7.6) on (7.5): 

u LM = (* — 1/6) — 1/12, u RtiJ = (i+ l/6)+l/12. 

(7.7) 

Then the solution by the piecewise-constant upwind 
scheme for this initial data after one time step is, 
by (7.4), 

U LJJ = UL ™ ~ ( 7 ' 8a ) 

and 

= RR,i,j C" x • (7.86) 

For the initial data (7.5), i.e., the data that are on 
the plane u = x, the exact solution is also given by 
(7.8). 

Thus, the bumpy initial data (7.7) is preserved by 
the piecewise-constant upwind scheme in the sense 
that the solution is given by first propagating the 
plane u — x exactly, and then superimposing the 
odd-even noise (7.6) back on. In fact, if we start 
with the data on the plane, the solution by the 
piecewise-constant upwind scheme turns it into a 
‘bumpy plane’ of type (7.7) after 8 iterations (for 
<j x = .2). Note that if we calculate the value at 
the center of each square by averaging the ‘bumpy’ 
values at the centroids of the two triangles, we get 
back the data on the plane, i.e., if the triangles are 


combined in pairs to form squares, then the cancel- 
lation effect again takes place. 

The above observation is consistent with the re- 
sult of Fourier analysis. To carry out this anal- 
ysis for the triangular mesh case, we must pair 
up the downward and upward pointing triangles 
so that the solution at each square looks the same 
as the solution at any other square. Consequently, 
each scheme has two amplification factors: princi- 
pal and spurious. The principal eigenfunctions of 
the (piecewise constant and piecewise linear) up- 
wind schemes turn out to be somewhat bumpy, 
while those for the centered and CE/SE schemes 
remain smooth. For the piecewise-constant upwind 
scheme, the spurious component is of order 0(h ) , 
but this error does not accumulate as we march to a 
fixed final time, and the piecewise constant upwind 
scheme retains its first-order accuracy. A similar re- 
mark also holds for the second-order upwind scheme 
(the spurious component is of order 0(h 2 )). 

To calculate the amplification factor, define (see 
Fig. 7.8) 


Then, the solution vector can be written as 

“Ij 1 = C -i,o u i-ij + c o,-i u «,j-i + c o,o u ;j 

+ C l,O u i+l,j CoyUj j +1 + - - - 

(7.10) 

Here each C; j is a 2 x 2 matrix. The two eigenval- 
ues of the matrix e~ IWx C_ 1 0 + e~ IWy C 0 x + . . . 
are the amplification factors of the corresponding 
scheme. The eigenvalue which approximates the ex- 
act amplification factor e 1 (- a * w *- (r yu>y) j s the prin- 
cipal amplification factor, and the other, the spuri- 
ous. 

Note that because the CE/SE scheme is odd-even 
decoupled, instead of pairing up the triangles, we 
only need to take two time steps. After calculating 
the eigenvalues (associated with two time steps), 
we need to take their square roots to obtain the 
amplification factor (per time step). 

The expressions for the amplification factors of 
the three schemes are very lengthy and are omitted 
here. The stability regions and accuracy compar- 
isons are shown below, however. 

Figures 7.11(a), (b), and (c) show the stability re- 
gions for the upwind, centered, and CE/SE schemes 
respectively. Here, the time step size for the upwind 
scheme is reduced considerably due to the stability 
limit along the two axes shown in Fig. 7.11(a). In 
practice, however, the maximum time step size is 
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not restricted too much compared with those of the 
centered and CE/SE schemes. 

For the next four plots, Figs. 7.12(a)-(c) and 
7.13, the wave numbers are w x = tt/16 and w y = 
ir/8. The four plots in Figs. 7.14(a)-(c) and 7.15 
are the same as the previous four except that 
w x — 7r/32 and w y = 7r/16. Observe that as in 
the rectangular-mesh case, the errors are of appro- 
priate order. 

Note that as shown by Figs. 7.13 and 7.15, the 
advantage in accuracy of the upwind scheme over 
the centered and CE/SE schemes is considerably 
less here compared to the rectangular-mesh case in 
Figs. 7.5 and 7.7. The upwind scheme also has the 
drawbacks of a somewhat large spurious component 
and a relatively small time step size (due to the 
flower-shape stability region). 

7.3. Accuracy comparison between rect- 
angular and triangular meshes. We can also 
compare the errors between the rectangular- and 
the triangular-mesh cases. First, observe that be- 
cause the triangular mesh is obtained by slicing each 
square into two triangles, the number of triangles 
are twice that of squares. To have the same num- 
ber of cells as the triangles, the squares must have 
a width of l/%/2- Thus, for the same number of 
cells, if the scheme is equally accurate between the 
square and the triangular mesh, then the phase er- 
rors in Figs. 7. 5-7.7 must be (\/2) 3 ~ 2.8 times 
the corresponding ones in Figs. 7.13-7.15, and the 
dissipation errors, (\/2) 4 = 4 times. 

For the upwind scheme, the maximum error in 
Fig. 7.7 is about 0.00004; that in Fig. 7.15 is also 
about 0.00004; i.e., there is no improvement in the 
triangular mesh case even though the number of 
cells is twice that of the quadrilateral mesh. Thus, 
for the same number of cells, the upwind scheme is 
more efficient on a rectangular mesh. 

Another drawback for the upwind scheme in the 
case of a triangular mesh is a rather large spuri- 
ous component. For a smooth data the spurious 
component is of order 0(h 2 ) whereas the error of 
the principal amplification factor per time step is 
0(h 3 ); here, h is the smallest edge length. To reach 
a fixed final time, we need a number of time steps 
proportional to 1 /h, and the error by the princi- 
pal amplification factor accumulates to 0(h 2 ). The 
spurious component is eventually damped out, and 
the corresponding error 0(h 2 ) does not accumulate. 
As a result, the (piecewise linear) upwind scheme 
retains its second-order accuracy for the case of a 
triangular mesh. 



(a) Upwind scheme, stability region, triangular 
mesh. 



(b) Centered scheme (i.e., extended N-T or coupled 
CE/SE). 



(c) CE/SE scheme. 

Fig. 7.11. Triangular-mesh stability regions. If the 
displacement vector based at the origin lies inside 
the stability region, the scheme is stable. The dark 
curve is produced by the principal amplification fac- 
tor; the light curve, the spurious one. The trian- 
gular cell and the reconstruction cell (hexagon) for 
the centered and CE/SE schemes are shown. 
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(a) Phase errors; triangular mesh 


(a) Phase errors; triangular mesh 


Dissipation Errors; flow along y=x 



Errors; flow along y=x 



Fig. 7.12. Errors for the triangular-mesh case; flow 
along y = x; w x — 7t/16 and w y = 7r/8. 


Dissipation Errors; flow along y=x 



(b) Dissipation errors 


Errors; flow along y=x 

0 . 00006 
0.00005 
0 . 00004 
0 . 00003 
0 . 00002 
0.00001 


0.1 0.2 0.3 0.4 

(c) Total errors 

Fig. 7.14. Errors for the triangular-mesh case; flow 
along y = x; w x = tt/ 32 and w y = 7r / 16 . 



Errors; average all directions 



Fig. 7.13. Total errors; average for all flow direc- 
tions; again w x = 7t/16 and w y = i r/8. 


Errors; average all directions 



Fig. 7.15. Total errors; average for all flow direc- 
tions; again w x = 7r/32 and w y = 7r/16. 
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For the centered and CE/SE schemes, the max- 
imum error in Fig. 7.7 is about 0.0002; that in 
Fig. 7.15 is about 0.00005; the improvement is by 
a factor of 4, more than the expected factor of 2.8. 
The maximum dissipation error in Fig. 7.7 is the 
value at a = CFL = 0, which is about 0.0001; that 
in Fig. 7.15 is about 0.000005, which is an improve- 
ment by a factor of 20, more than the expected 
factor of 4. Thus, between a quadrilateral and a 
triangular mesh with the same number of cells, the 
centered and CE/SE schemes are more efficient on 
the triangular mesh. 

Notice that the above analysis is linear. It pro- 
vides useful information on stability and accuracy. 
However, the analysis no longer holds when there is 
a solid wall or when a limiter function is employed. 
In such cases, we resort to numerical experiments. 

8. Numerical results. Results for the ID as 
well as the 2D quadrilateral and triangular meshes 
are shown below. 

8.1 Results for the ID case. In the following 
ID numerical examples, 7 = 1.4 and the CFL num- 
ber is 0.9. Unless otherwise stated, the number of 
mesh points is 200. 

The first numerical test, used by Sod (1978), is 
the Riemann problem 

(pL i H L i Pl') ( 1 , 0 , 1 ), 

( PR, U R , PR ) = (0.125,0,0.1). 

The final time is t = 0.2. The solid line repre- 
sents the exact solution. The results are shown in 
Fig. 8.1. 

The second problem, due to Shu and Osher 
(1989), has several extrema in the smooth regions. 
In the interval — 5 < x < 5, a moving Mach 3 shock 
interacts with sine waves in density as described by 
the following initial conditions: 

, x _ / (3.857, 2.629, 10.333) if x < -4, 

(P, u,p) - | ( 1 + 0 2 s i n 5a;, 0, 1) otherwise. 

The final time is t — 1.8. Figure 8.2 shows the re- 
sults with 800 mesh points. Since the exact solution 
is not known, the solid line represents the solution 
for 1600 cells via a uniformly second-order accurate 
upwind scheme described in Huynh (1995a). 

Note that for these two problems, the upwind 
results are more accurate than the N-T and CE/SE 
results. Between the N-T and CE/SE solutions, the 
formers are very slightly more dissipative. 



0.0 0.2 0.4 0.6 0.8 1.0 


X 

Fig. 8.1 Sod’s problem. 



- 3 - 2-10 1 2 3 


x 

Fig. 8.2 Shu and Osher’s problem (800 mesh points). 



0.4 0.6 0.8 1.0 


X 

Fig. 8.3 Slow-moving shock problem. Final time = 2.2 
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density 



Fig. 8.4 Straightforward nonreflecting boundary 
condition works for the CE/SE as well as the up- 
wind and N-T schemes. 



Fig. 8.5. Oblique shock; 80 x 20 rectangular and 
40 x 10 x 4 triangular meshes. 


The next test is the slow-moving shock problem 
for which the upwind scheme generates oscillations. 
The initial condition is a velocity of 0.05 superim- 
posed on a Mach 3 normal shock: 

(Pl,ul,Pl) = (0.24,2.2238,0.09), 

(Pr,ur, P r ) = (0.9256,0.6137,0.9299). 

Figure 8.3 shows the results for the final time t = 3. 
Here, the upwind solution oscillates. Note that the 
entropy condition (4.49) employed helps reduce os- 
cillations considerably compared with a typical im- 
plementation of Roe’s scheme or (4.48). The N-T 
and CE/SE solutions have no oscillations. An up- 
wind step that generates no oscillations for this 
problem can be found in Wada and Liou (1997). 

The last ID test concerns nonreflecting bound- 
ary conditions, which let waves leave the domain 
of computation. The boundary condition employed 
below is an obvious one: the value at the ghost cell 
is set to be that of the immediate (interior) neigh- 
bor, and the slope is set to zero (see also Chang et 
al. 1999). The results for Sod’s problem at time 
.28, .30, .40, and .60 are shown in Fig. 8.4. Here, 
since the CE/SE solutions are again essentially the 
same as the N-T solutions, they are omitted. The 
nonreflecting boundary condition works well for all 
schemes: there are no reflecting waves. 


8.2 Results for the 2D case. As explained at 
the beginning of §7, the extended N-T scheme us- 
ing the CE/SE approach to extension is called the 
centered scheme here — as opposed to the upwind 
scheme. If the slopes are carried along, the corre- 
sponding extension is called the CE/SE scheme. 

For 2D test problems, the solutions are obtained 
by the author’s codes except in the case of the 
CE/SE scheme on a triangular mesh where the 
CE/SE code, kindly provided by Dr. Xiao- Yen 
Wang, is employed. Because a unified CFL number 
for both the triangular and quadrilateral meshes is 
not readily available, we provide information on the 
time steps, which are generally chosen to be close 
to the largest time step possible. 

The schemes discussed here are designed to solve 
unsteady problems. The 2D test cases chosen be- 
low, however, are steady-state problems because the 
exact solutions are either known or have certain fea- 
tures which are useful for the purpose of compari- 
son. 

The first 2D test is the oblique shock problem 
on the domain [0,4] x [0,1] (Yee et al. 1983). The 
condition at the inflow boundary is (p, M,i>,p) = 
(1., 2.9,0., 1./1.4) and that at the top boundary, 
(1.7, 2.6193, —.5063, 1.5282). The rectangular mesh 
is 80 x 20. The triangular mesh is obtained by 
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slicing each rectangle (along the two diagonals) 
of a 40 x 10 rectangular mesh into four trian- 
gles. Thus, each mesh has 1600 cells. The solu- 
tions are shown in Fig. 8.5; here, pressure at cell 
centers at a fixed y location are plotted; for the 
rectangular-mesh case y = .475; for the triangular- 
mesh case y = .45. Note that all schemes smear the 
oblique shocks. The upwind solution is slightly less 
smeared than the centered and CE/SE solutions for 
the rectangular-mesh case, but it oscillates for the 
triangular-mesh case, while the other two solutions 
do not. Also note that again, the solutions by the 
centered and CE/SE schemes are nearly identical. 
For these two schemes, the triangular-mesh solu- 
tions are also nearly the same as the rectangular- 
mesh solutions. 

Notice that the CE/SE solutions shown in 
Fig. 8.5 here are considerably different from the so- 
lutions by Chang et al. (1999) and Zhang et al. 
(2002), where the shocks are resolved by only one 
point. The reason is that the mesh is different. The 
result for the 80 x 20 rectangular mesh here has 
also been verified by Dr. Xiao- Yen Wang using her 
CE/SE code (private communication). 

Next, we discuss the time step sizes and con- 
vergence. In the quadrilateral-mesh case, for the 
upwind scheme, with At = 0.0106, the solution 
converges to machine accuracy. For the centered 
scheme, with At = 0.0077, the solution converges 
to machine accuracy. For the CE/SE scheme, the 
solution does not converge well for most At, e.g., 
with At = 0.0064, the maximum residual goes down 
three orders of magnitude, and with At = 0.005, 
four orders of magnitude (the other two schemes 
converge to machine accuracy for these time step 
sizes). It appears that because the scheme is odd- 
even decoupled, the solution does not converge well. 
Note that for the centered and CE/SE schemes, 
the residual is calculated using the data at even 
time levels only (as discussed at the end of §1.3). 
In the triangular-mesh case, the solution converges 
to machine accuracy for the upwind scheme with 
At = 0.0078, and the centered scheme, At = 0.008. 
For the CE/SE scheme, with At = 0.0074, the max- 
imum residual goes down three orders of magnitude. 

Figure 8.6(a), (b), and (c) show the solutions on 
an unstructured triangular mesh. (The mesh for 
this case as well as that in Fig. 8.9 were kindly 
provided by Dr. Philip Jorgenson.) Here, the up- 
wind solution is slightly less smeared compared with 
the other two. Also note that this upwind solution 
(Fig. 8.6(a)) improves considerably compared with 
the oscillatory solution in Fig. 8.5; the minimum 


pressure here is .698 whereas that for the 40 x 10 x 4 
triangular mesh, .612 (and that for the y — .45 slice 
in Fig. 8.5, .653); the exact minimum pressure is 
.714. 

Since solutions by the CE/SE schemes are essen- 
tially the same as those by the centered schemes, 
they are omitted from here on. 

The next problem is the nozzle problem (Ver- 
hoff 1985). Here, the domain is [0,3] x [0, 1] except 
that for x in [1,2], the bottom wall is defined by 
y = .lsin 2 (ar — 1 )tt. The flow is subsonic, and the 
boundary conditions are standard. At inflow, we 
set the total pressure po = 1; stagnation speed of 
sound ao = 1; and Vq = 0: pressure, however, is ex- 
trapolated from the interior. At outflow, pressure 
is set to .843, which results in an inflow Mach num- 
ber of .5; the other three variables are extrapolated 
from the interior. Since there are no shocks in the 
solution, the average slope is employed. 

The 64 x 16 quadrilateral mesh and the solutions 
are shown in Fig. 8.7(a) and (b) respectively. Here, 
in Fig. 8.7(b), the line is the solution by the cen- 
tered scheme on a 128 x 32 mesh. Observe that the 
solution by the upwind scheme is slightly more ac- 
curate: the peak is higher, and the valleys are lower 
than those of the centered scheme. 

For the same problem, the 32 x 8 x 4 triangular 
mesh and the solutions are shown in Fig. 8.8(a) and 
(b) respectively. Here, in Fig. 8.8(b), the solution 
by centered scheme is slightly more accurate and 
maintains symmetry better than that of the upwind 
scheme. Also note the accuracy at the lower wall 
near the outflow boundary. 

The last test is the transonic flow over a circu- 
lar bump. The boundary conditions are the same 
as those of the above nozzle problem except that 
the outflow pressure is .736. The mesh has 2466 
cells. For the upwind scheme, the time step is .011. 
After 5000 iterations, the residual decreases from 
.938E-1 to .918E-4, but after another 24000 itera- 
tions, the residual is still .532E-4, i.e., the solution 
does not converge well. A time step of 0.008 does 
not improve convergence. For the centered scheme, 
the time step is .013 (.014 does not converge). After 
10000 iterations, the residual decreases from .157E0 
to .822E-8, and after another 10000 iterations, the 
residual reaches machine error of E-14. 

9. Conclusions and discussion. The numer- 
ical tests here are clearly preliminary; more tests 
need to be performed. 

In summary, the second-order accurate upwind, 
the Lax-Friedrichs type second-order accurate cen- 
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(a) Upwind scheme 
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(b) Centered scheme 
(c) CE/SE scheme 

Fig. 8.6. Pressure contours, oblique shock. 


Fig. 8.8(a) Triangular mesh (32 x 8 x 4); nozzle 
problem. 
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Fig. 8.7(a) Quad mesh (64 x 16); nozzle problem,. Fl & 8 ‘ 8 ( b ) Triangular-mesh results; nozzle problem. 


(a) Upwind scheme 


0.0 0.5 1.0 1.5 2.0 2.5 3.0 


(b) Centered scheme 

Fig. 8.9. Mach contours, circular bump on a wall. 


Fig. 8.7(b) Quad-mesh results; nozzle problem. 
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tered (or the extended N-T scheme using the 
CE/SE approach to extension), and the CE/SE 
(e = 1/2) schemes were presented in a framework 
that facilitates their comparison. These schemes all 
use piecewise-linear reconstructions (MUSCL inter- 
polants). The key difference is that the centered 
and CE/SE schemes avoid the upwind step by em- 
ploying reconstruction cells that are roughly twice 
as big as the original cells. 

The schemes employed have several properties 
in common. They are all of finite-volume type, 
conserve fluxes locally and globally, are equally di- 
mensional splitting (or non-splitting), can capture 
shocks, are only first-order accurate near extrema 
when the slope is defined by a weighted average 
or a limiter function, and can handle nonreflecting 
boundary conditions. 

Fourier stability and accuracy analyses were car- 
ried out for these schemes for the standard ID and 
the 2D quadrilateral mesh cases. In the nonstan- 
dard case of a triangular mesh, the triangles were 
paired up for the analyses of the upwind and the 
centered schemes. Among the three schemes, on a 
quadrilateral mesh, the upwind scheme is more ac- 
curate. On a triangular mesh, however, accuracy 
among the three schemes appears to be compara- 
ble. Comparing the same scheme on the two meshes 
with the same number of cells, the upwind scheme is 
more efficient on a quadrilateral mesh, and the cen- 
tered and CE/SE schemes, on a triangular mesh. 

The centered scheme has numerous advantages 
over the CE/SE scheme: it produces essentially the 
same solution, runs faster, converges better, and 
requires only one third of the storage for flow vari- 
ables. The capability of adjusting dissipation of the 
CE/SE schemes, however, has not been fully ex- 
plored. 

Compared with the upwind scheme, the centered 
scheme is more stable, conceptually simpler (the ge- 
ometry involved is slightly more complicated, how- 
ever) , and requires less storage (for an unstructured 
mesh). Besides providing an alternative to the up- 
wind choice, the centered scheme can be very use- 
ful when the equations become complicated, and 
the upwind model either has not been derived or is 
no longer appropriate (the equations are no longer 
hyperbolic). 

The development, analysis, and comparison of 
these schemes with viscous terms remain to be ex- 
plored. (A discussion of the CE/SE scheme for vis- 
cous flows can be found in Chang et al. (1995).) 
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